# User Manual Movgrid

User Manual: Pdf

Open the PDF directly: View PDF .

Page Count: 11

User’s Manual for the Movgrid Toolbox

The MatMOL Group (2009)

Introduction

The Movgrid toolbox consists of a set of functions included in MatMOL which allows the

user to solve systems described by partial differential equations (PDEs) using the Dynamic

Regridding approach. This dynamic regridding is carried out based in the equidistribution

principle. The algorithm for the moving grid was originally proposed in [1, 3]. A description

(including some benchmark examples) of the algorithm employed in the MatMOL toolbox is

given in [2].

The Movgrid Toolbox by the Example

In this section we will show how to use the Movgrid toolbox to solve the following hypothetic

system of PDEs:

∂u1

∂t =−∂f1(u3)

∂z +ε∂2u1

∂z2+k1u2(1)

∂u2

∂t =−∂f2(u2, u3)

∂z +ε∂2u2

∂z2(2)

∂u3

∂t =−∂f3(u1, u2, u3)

∂z +ε∂2u1

∂z2(3)

where

f1=u3;f2= (γ−1)u3−(γ−3)u2

2

2u3

;f3=γu3−(γ−1)u2

2

2u1u2

u1

with boundary and initial conditions of the form:

∂u1(0, t)

∂z = 0; ∂u1(1, t)

∂z = 0; (4)

u2(0, t) = 0; u2(1, t) = 0; (5)

∂u3(0, t)

∂z = 0; ∂u3(1, t)

∂z = 0; (6)

u1(z, 0) = 1if 0≤z≤0.5

0.125 if 0.5< z ≤1;u2(z, 0) = 0;

u3(z, 0) = 2.5if 0≤z≤0.5

0.25 if 0.5< z ≤1(7)

The structure of the software of the Movgrid in Matlab c

is as follows:

1

•A set of Matlab c

functions that performs automatic operations for computing the differ-

ent terms in the moving grid algorithm (see [2] for details). It should be mention that this

functions are independent of the problem and they are included in the MatMOL toolbox.

•A main program that allows the user to specify the space and time domains, the PDE

parameters, the initial conditions, the type of discretization scheme, the type of monitor

function, the tuning parameters of the adaptive grid, etc. This ﬁle is problem dependent.

In the sequel we will explain how to construct the different elements of this ﬁle which

we will call Main_file.m.

•A set of MATLAB functions where a user input is needed for deﬁning

a) The right-hand side of the PDEs (1)-(3) is included in right_vect.m

b) The boundary conditions (4)-(6) are included in Bcintroduct.m

c) Possibly the ﬂux function if the PDE is in conservation form

d) Possibly the Jacobian of the ﬂux function if the PDE is in conservation form

These ﬁles are also problem dependent

Let us start the description of the implementation with the Main_file.m.

The structure of Main file.m

First, the user should indicate the global variables that will be employed by the other ﬁles in

the toolbox

% Global variables

global zL zR n ne X eq choice

global alpha kappa mu tau

global u z

global gam eps k1

The ﬁrsts three lines of global variables are common to all the problems. The last one is

for the PDE parameters that the user wants to pass to the right_vect.m function.

No the user can deﬁne, for instance, the parameters which are common to other ﬁnite

difference schemes like the temporal conditions, the parameters of the PDEs and the spatial

grid

% Temporal conditions

t = [0 .15 .23 .28 .32];

% PDE parameters

ne = 3; % Number of PDE equations

eps = 10ˆ(-3);

gam = 1.4;

k1 = 1;

2

% Initial spatial grid z : if n is the number of moving

% grid points, then dim(z) = n+2

zL = 0; % Initial point of the grid

zM = 0.5; % A middle point of the grid

zR = 1; % Final point of the grid

nLM = 21; % Number of points in the first interval

nMR = 61; % Number of points in the second interval

dzLM = (zM-zL)/(nLM-1);

dzMR = (zM-zR)/(nMR-1);

z = [zL:dzLM:zLM zLM+dzMR:dzMR:zMR]’; % The grid

n = length(z)-2; % Moving points

where ne is the number of partial differential equations. A middle point (zM) was included to

introduce the initial conditions.

After this, the user can deﬁne the parameters for the moving grid algorithm (see [2] for

details).

% Moving grid parameters

alpha = 0.05;

kappa = 1;

mu = kappa*(kappa+1);

tau = 1e-4;

In the next block of the m-ﬁle, the initial conditions (7) will be implemented:

% Initial conditions

for i=1:n+2,

if (z(i) <= zM)

u(i,1) = 1;

u(i,2) = 0;

u(i,3) = 2.5;

else

u(i,1) = 0.125;

u(i,2) = 0;

u(i,3) = 0.25;

end

end

Now the user has to specify the variables that function JPat will employ to implement the

sparsity pattern of the Jacobian in X. A tensor eq of dimension ne×3×ne, with nebeing the

number of PDEs, is created based on the PDEs (1)-(3) and the chosen ﬁnite difference scheme.

Consider a three point centered ﬁnite difference scheme for the ﬁrst derivative and a ﬁve

point centered ﬁnite difference scheme for the second derivative. Let us now show how the

tensor eq is constructed in Matlab c

% Global variables for JPat

X = zeros(n*(ne+1),n*(ne+1));

% Tensor eq

3

eq(:,:,1) = [-2 2 1 ; 0 0 1 ; -1 1 1];

eq(:,:,2) = [ 0 0 0 ; -2 2 1 ; -1 1 1];

eq(:,:,3) = [-1 1 1 ; -1 1 1 ; -2 2 1];

In the construction of eq one can note the following:

•Matrix eq(:,:,k) describes the pattern of points used in equation kin order to take

into account the chosen ﬁnite difference schemes.

•The line jin matrix eq(:,:,k) refers to the presence of variable j(uj) in equation k.

•In order to ﬁll the elements of the matrix eq(:,:,k) we look at the ﬁnite difference

schemes.

–For the computation of the second derivative in the point i, the following points

will be involved: i−2,i−1,i,i+ 1 and i+2(ﬁve point centered scheme). The

extremes points in this case are i−2and i+2.

–For the computation of the ﬁrst derivative in the point i, the following points will

be involved: i−1,iand i+1(three point centered scheme). The extremes points

in this case are i−1and i+1

The ﬁrst column of matrix eq(:,:,k) refer to the left extreme of the ﬁnite difference

scheme (in the case of the second derivative i−2so −2, in the case of the ﬁrst derivative

i−1so −1). The second column of matrix eq(:,:,k) refer to the right extreme of

the ﬁnite difference scheme (in the case of the second derivative i+ 2 so 2, in the case of

the ﬁrst derivative i+ 1 so 1). The third column indicates either if the variable appears

implicity in the equation (in which case it takes the value 1) or not (in which case it takes

the value 0).

•If we look at the ﬁrst equation -Eqn (1)-, we see that:

–The ﬁrst variable (u1) is present only via the second derivative ∂2u1

∂z2term. Since

for the second derivative we have chosen a ﬁve point centered scheme with ex-

tremes i−2and i+2, we have:

eq(1,:,1) = [-2 2 1]

–The second variable (u2) appears only without derivative (k1u2) so :

eq(2,:,1) = [0 0 1]

–The third variable (u3) appears only in the ﬁrst derivative term through ∂f1(u3)

∂z .

Since for the ﬁrst derivative we have chosen a three point centered scheme with

extremes i−1and i+1, we have:

eq(3,:,1) = [-1 1 1]

•If we look at the second equation -Eqn (2)-, we see that:

4

–The ﬁrst variable (u1) does not appear so:

eq(1,:,2) = [0 0 0]

–The second variable (u2) appears in the second derivative term and in the ﬁrst

derivative term through function f2(u2, u3). In this case we have to choose the

set that include all the points, so:

eq(2,:,2) = [-2 2 1]

Since [−1,1] is included in [−2,2].

–The third variable (u3) appears only in the ﬁrst derivative term through f2(u2, u3)

so:

eq(3,:,2) = [-1 1 1]

•Similarly we can construct eq(:,:,3) by looking at the third equation -Eqn (3)-.

Example 1 (Another example for constructing tensor eq(:,:,k))Consider the system:

∂u1

∂t =D∂2u1

∂z2−v∂u1

∂z +u2

∂u2

∂t =−v∂u2

∂z

Compute the tensor eq(:,:,:) if for the ﬁrst derivative a ﬁve point biased upwind ﬁnite

difference scheme (i−3,i−2,i−1,iand i+1) is employed, while for the second derivative

a ﬁve point centered ﬁnite difference (i−2,i−1,i,i+1and i+2) is preferred.

If we look at the ﬁrst equation, we see that:

•The ﬁrst variable (u1) is present via the ﬁrst and second derivative terms. In this case

we have to choose the set that include all the points, so:

eq(1,:,1) = [-3 2 1]

Since the set [−3,2] includes both the set [−3,1] from the ﬁve point biased upwind and

the set [−2,2] from ﬁve point centered.

•The second variable (u2) appears only without derivative (u2) so :

eq(2,:,1) = [0 0 1]

If we look at the second equation, we see that:

•The ﬁrst variable (u1) does not appear so:

eq(1,:,2) = [0 0 0]

•The second variable (u2) is present only via the ﬁrst derivative term so:

eq(2,:,2) = [-3 1 1]

5

In the next block of Main file.m to solve system (1)-(7) the user will choose the monitor

function employed by the adaptive grid. Three different possibilities are available:

•’der1’. Monitoring based on the ﬁrst derivative of the ﬁeld.

m(i) = v

u

u

tα+1

ne

ne

X

j=1

(u(i+ 1, j)−u(i, j))2

(z(i+ 1) −z(i))2

•’der2’. Monitoring based on the second derivative of the ﬁeld.

m(i) = v

u

u

tα+1

ne

ne

X

j=1

∂2u(i+1,j)

∂z2+∂2u(i,j)

∂z2

2

•’derf1’. Monitoring based on the ﬁrst derivative of the ﬂux.

m(i) = v

u

u

tα+1

ne

ne

X

j=1

(fl[u(i+ 1, j)] −f l[u(i, j)])2

(z(i+ 1) −z(i))2

This choice is easily implemented in Matlab c

, for instance:

% Monitoring choice

choice = ’der1’;

Usually when we are considering a conservation equation (like in this example) ’der1’

and ’derf1’ work better than ’der2’. On the contrary when trying to solve wave equations,

’der2’ will be, usually, preferable.

It should be mentioned here that, in some problems, it may be useful to perform an adapta-

tion of the initial grid. Some examples indicating how to do this were included in the MatMOL

toolbox (see, for instance, fischkol main.m or burger main.m ﬁles).

In the moving grid algorithm, the state variables and the coordinates of the internal (mov-

ing) grid points are taken together to solve the problem. The new vector xconstructed from

the state variables and grid coordinates is of the form:

x= [u1

1, u2

1, ..., une

1, z1, u1

2, u2

2, ..., une

2, z2, ..., u1

n, u2

n, ..., une

n, zn]T

where the subindices indicate the position in the grid and the superindices denote the ﬁeld

considered. In this way, uj

irefers to the value of the ﬁeld ujin the grid point zi.

The next block in the script Main file.m will assemble state variables and grid coordi-

nates in x:

% Global initial condition

for jj = 1:ne,

x(jj:ne+1:n*(ne+1)) = u(2:n+1,jj);

end

x(ne+1:ne+1:n*(ne+1)) = z(2:n+1);

6

Finally we choose the parameters of the ODE solver and call it using, for instance, the

following instructions

% ODE Solver parameters

options = odeset(’RelTol’,1e-6,’AbsTol’,1e-6,’Mass’,@mass,...

’MStateDependence’,’strong’,’JPattern’,JPat,...

’MvPattern’,MvPat,’MassSingular’,’no’);

% Call to the ODE solver

[tout, yout] = ode15s(@right_vect,t,x,options);

After this, the user can include a block for visualizing the results.

The structure of right vect.m

After applying the moving grid approach to eqns (1)-(2) a semidiscrete system of ODEs is

obtained. The general form of this system is [2]:

A(x, t)dx

dt =b(x, t)

where A(x, t)is the mass matrix and b(x, t)is the discrete version of the right hand side terms

of eqns (1)-(2) with xbeing the global coordinates as indicated in Main file.m. This is, a

vector which includes the state variables and the grid coordinates.

In this ﬁle the user will deﬁne the right hand side terms (b(x, t)). These terms are imple-

mented in matrix udot(:,k) where kindicates the number of the equation. In this ﬁle, the

user will also choose the ﬁnite difference schemes for the computation of the spatial derivatives.

It should be noted that vector xcontains only the interior (moving) grid points, nevertheless,

b(x, t)also makes use of the boundary points z0,zn+1,uj

0,uj

n+1 with j= 1, .., ne. The

structure of b(x, t)∈R(ne+1)n×1is:

b(x, t) = du1

1

dt ,du2

1

dt , ..., dune

1

dt , g1,du1

2

dt ,du2

2

dt , ..., dune

2

dt , g2, ..., du1

n

dt ,du2

n

dt , ..., dune

n

dt , gnT

where

gi=1

Mi−µ

∆zi+1

+1+2µ

∆zi

−µ

∆zi−1−1

Mi−1−µ

∆zi

+1+2µ

∆zi−1

−µ

∆zi−2

For details see [1, 3].

In the right vect.m function the user will ﬁrst declare the global variables employed

also in other functions of the toolbox:

% Global variables

global n ne choice

global dltz mu

global gam eps k1

7

The second block will make a call to the function which implements the boundary condi-

tions (this function, as we will show in the next section, will also separate the state variables

and the node positions).

% Separate state variables and node positions and implement the BCs

[u z] = Bcintroduct(t,x);

In the next block the user will compute the monitoring function (by calling the Matlab c

function monitor.m) and the right vector g(by calling function rightmon.m)

% Compute the monitoring function mon(i)

mon = monitor(u,z,choice,t);

% Compute the right vector g(n) of the moving grid equation

g = rightmon(mon,dltz,mu);

The blocks presented so far are independent of the problem (except the last line of global

variables).

After this, the user will deﬁne the nonlinear functions, choose the ﬁnite difference scheme

and compute the spatial derivatives of the ﬁelds and nonlinear functions:

% Nonlinear functions

f1 = u(:,3);

f2 = (gam-1)*u(:,3)- ((gam-3)*u(:,2).ˆ2)./(2*u(:,3));

f3 = (gam*u(:,3) - (gam-1)*u(:,2).ˆ2./(2*u(:,1))).*u(:,2)./u(:,1);

% Compute the finite difference matrices

D1 = matfd(z, ’non_uni’, ’1st’, 3, ’centered’);

D2 = matfd(z, ’non_uni’, ’2nd’, 5, ’centered’);

% Field second derivatives

u1zz = D2*u(:,1);

u2zz = D2*u(:,2);

u3zz = D2*u(:,3);

% Nonlinear functions first derivatives

f1z = D1*f1;

f2z = D1*f2;

f3z = D1*f3;

The implementation of the right hand side terms is carried out in the next block:

% Temporal derivatives (RHS terms)

udot(:,1) = -f1z + eps*u1zz + k1*u(:,2);

udot(:,2) = -f2z + eps*u2zz;

udot(:,3) = -f3z + eps*u3zz;

Finally, the last block, which is independent of the problem, is for the assembly of the

global right vector b(x, t)

8

% Assemble the global right vector

for jj = 1:ne

xt(jj:ne+1:n*(ne+1),1) = udot(2:n+1,jj);

end

xt(ne+1:ne+1:n*(ne+1),1) = g;

The structure of Bcintroduct.m

In this function the user must separate the state variables from the grid positions and implement

the boundary conditions. The boundary conditions need to be speciﬁed in Dirichlet form. This

means that for Neumann and Robin type a transformation have to be performed.

The ﬁrst part (separate state variables from grid positions) is the same for all problems and

it is easily carried out with the following Matlab c

lines

% Separate dependent variables and node positions

for jj = 1:ne

u(2:n+1,jj) = x(jj:ne+1:n*(ne+1));

end

z = [zL x(ne+1:ne+1:n*(ne+1))’ zR];

For the implementation of the boundary conditions in Dirichlet form, two possibilities are

available:

•To compute the value of the ﬁeld in the boundary by hand. If, for instance, we have the

following boundary condition at the right boundary:

a1

∂uk(zR)

∂z +a0uk(zR) = f(t)(8)

we can approximate the ﬁrst derivative by a simple ﬁnite difference scheme so that

∂uk(zR)

∂z =u(n+ 2, k)−u(n+ 1, k)

z(n+ 2) −z(n+ 1)

Substituting this expression into equation (8) and after some algebraic manipulations we

obtain:

u(n+ 2, k) =

f(t) + a1

u(n+ 1, k)

z(n+ 2) −z(n+ 1)

a0+a1

z(n+ 2) −z(n+ 1)

•To use the MatMOL function BC_dir. This function is called as follows:

u = BC_dir(bc_form, u, uin, D1, pos, v, Dif)

% Input parameters:

% bc_form : Type of the boundary condition ’dir’ for Dirichlet

% boundary type, ’neu’ for Neumann type and ’rob’

% for Robin type.

% u : Value of the field in the spatial domain

% uin : Value of the non homogeneous part of the boundary.

9

% For homogeneous boundary conditions, set this

% value to 0.

% D1 : Finite difference matrix for the first derivative

% pos : Position of the boundary condition: ’begin’ if the

% boundary condition refers to the first point of

% the spatial domain, ’end’ if it refers to the last

% point.

% v : Fluid velocity. This parameter is only necessary

% in the case of Robin boundary conditions,

% otherwise it can be set to 0

% Dif : Difussivity. This parameter is only necessary in

% the case of Robin boundary conditions, otherwise

% it can be set to 0

%

% Output parameters:

% ubc : Value of the field in the boundary (Dirichlet

% form)

Using this function, boundary condition (8) can be computed as

Dz = matfd (z, ’non_uni’, ’1st’, 3, ’upwind’,v);

u(n+2,k) = BC_dir(’rob’, u(:,k), f/a0, Dz, ’end’, a1, a0);

Boundary conditions (4)-(6) are implemented in Matlab c

as follows

% Implement the BCs

% Computation of the derivative matrix

Dz = matfd (z, ’non_uni’, ’1st’, 5, ’centered’);

% Boundary conditions for u1

u(1,1) = BC_dir(’neu’,u(:,1),0,Dz,’begin’,0,0);

u(n+2,1) = BC_dir(’neu’,u(:,1),0,Dz,’end’,0,0);

% Boundary conditions for u2

u(1,2) = 0;

u(n+2,2) = 0;

% Boundary conditions for u3

u(1,3) = BC_dir(’neu’,u(:,3),0,Dz,’begin’,0,0);

u(n+2,3) = BC_dir(’neu’,u(:,3),0,Dz,’end’,0,0);

Note that other ﬁnite difference schemes can be used (when choosing D1) to approximate

the boundary conditions.

Remark: The boundary conditions must be always introduced, even in the cases where

mathematically it is not necessary. Consider, for instance, an equation which is purely convec-

tive. From the mathematical point of view, only one boundary condition is required to solve

it. Nevertheless, with the moving grid algorithm, both boundary conditions (at the beginning

and at the end) have to be speciﬁed. It is necessary, then, to intoduce a “ﬁctitious” boundary

condition. From the point of view of time integration, the most ﬂexible way of deﬁning such

10

conditions is to make the ﬁrst spatial derivative of the corresponding function to be constant. If

we consider the ﬁrst point of the spatial domain, this translates into:

∂u(1, k)

∂z =∂u(2, k)

∂z

References

[1] J. G. Blom and P. A. Zegeling. Algorithm 731: A moving-grid interface for systems

of one-dimensional time-dependent partial differential equations. ACM Transactions on

Mathematical Software, 20:194–214, 1994.

[2] P. Saucez, L. Some, and A. Vande Wouwer. Matlab implementation of a moving grid

method based on the equidistribution principle. Applied Mathematics and Computation,

d.o.i. 10.1016/j.amc.2009.07.34, 2009.

[3] J. G. Verwer, J. G. Blom, R. M. Furzeland, and P. A. Zegeling. A moving-grid method

for one-dimensional pdes based on the method of lines. In J.E. Flaherty, P.J. Paslow, M.S.

Shephard, and J.D. Vasilakis, editors, Adaptive Methods for Partial Differential Equations,

pages 160–175. SIAM, Philadelphia, 1989.

11