User Manual Movgrid
User Manual: Pdf
Open the PDF directly: View PDF .
Page Count: 11
Download | |
Open PDF In Browser | View PDF |
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 on 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 ∂f1 (u3 ) ∂ 2 u1 =− + ε 2 + k1 u2 ∂t ∂z ∂z ∂f2 (u2 , u3 ) ∂ 2 u2 ∂u2 =− +ε 2 ∂t ∂z ∂z ∂f3 (u1 , u2 , u3 ) ∂ 2 u1 ∂u3 =− +ε 2 ∂t ∂z ∂z (1) (2) (3) where f1 = u3 ; (γ − 3)u22 f2 = (γ − 1)u3 − ; 2u3 (γ − 1)u22 u2 f3 = γu3 − 2u1 u1 with boundary and initial conditions of the form: ∂u1 (0, t) = 0; ∂z u2 (0, t) = 0; ∂u1 (1, t) = 0; ∂z (4) u2 (1, t) = 0; (5) ∂u3 (0, t) ∂u3 (1, t) = 0; = 0; ∂z ∂z 1 if 0 ≤ z ≤ 0.5 u1 (z, 0) = ; u2 (z, 0) = 0; 0.125 if 0.5 < z ≤ 1 2.5 if 0 ≤ z ≤ 0.5 u3 (z, 0) = 0.25 if 0.5 < z ≤ 1 The structure of the software of the Movgrid in Matlab c is as follows: 1 (6) (7) • A set of Matlab c functions that performs automatic operations for computing the different 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 file is problem dependent. In the sequel we will explain how to construct the different elements of this file which we will call Main_file.m. • A set of MATLAB functions where a user input is needed for defining 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 flux function if the PDE is in conservation form d) Possibly the Jacobian of the flux function if the PDE is in conservation form These files 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 files 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 firsts 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. Now the user can define, for instance, the parameters which are common to other finite difference schemes like the temporal conditions, the parameters of the PDEs and the spatial grid % Temporal conditions t = [0 .15 .23 .28 .32]; % PDE ne = eps = gam = k1 = parameters 3; 10ˆ(-3); 1.4; 1; % Number of PDE equations 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 define 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-file, 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 matrix eq of dimension ne × 3 × ne , with ne being the number of PDEs, is created based on the PDEs (1)-(3) and the chosen finite difference scheme. Consider a three point centered finite difference scheme for the first derivative and a five point centered finite difference scheme for the second derivative. Let us now show how the matrix eq is constructed in Matlab c % Global variables for JPat X = zeros(n*(ne+1),n*(ne+1)); % Matrix 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 k in order to take into account the chosen finite difference schemes. • The line j in matrix eq(:,:,k) refers to the presence of variable j (uj ) in equation k. • In order to fill the elements of the matrix eq(:,:,k) we look at the finite 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 (five point centered scheme). The extremes points in this case are i − 2 and i + 2. – For the computation of the first derivative in the point i, the following points will be involved: i − 1, i and i + 1 (three point centered scheme). The extremes points in this case are i − 1 and i + 1 The first column of matrix eq(:,:,k) refer to the left extreme of the finite difference scheme (in the case of the second derivative i − 2 so −2, in the case of the first derivative i − 1 so −1). The second column of matrix eq(:,:,k) refer to the right extreme of the finite difference scheme (in the case of the second derivative i + 2 so 2, in the case of the first 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 first equation -Eqn (1)-, we see that: – The first variable (u1 ) is present only via the second derivative ∂ 2 u1 ∂z 2 term. Since for the second derivative we have chosen a five point centered scheme with extremes i − 2 and i + 2, we have: eq(1,:,1) = [-2 2 1] – The second variable (u2 ) appears only without derivative (k1 u2 ) so : eq(2,:,1) = [0 0 1] – The third variable (u3 ) appears only in the first derivative term through ∂f1 (u3 ) ∂z . Since for the first derivative we have chosen a three point centered scheme with extremes i − 1 and i + 1, we have: eq(3,:,1) = [-1 1 1] • If we look at the second equation -Eqn (2)-, we see that: 4 – The first 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 first 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 first 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 matrix eq(:,:,k)) Consider the system: ∂ 2 u1 ∂u1 ∂u1 =D 2 −v + u2 ∂t ∂z ∂z ∂u2 ∂u2 = −v ∂t ∂z Compute the matrix eq(:,:,:) if for the first derivative a five point biased upwind finite difference scheme (i − 3, i − 2, i − 1, i and i + 1) is employed, while for the second derivative a five point centered finite difference (i − 2, i − 1, i, i + 1 and i + 2) is preferred. If we look at the first equation, we see that: • The first variable (u1 ) is present via the first 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 five point biased upwind and the set [−2, 2] from five 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 first variable (u1 ) does not appear so: eq(1,:,2) = [0 0 0] • The second variable (u2 ) is present only via the first derivative term so: eq(2,:,2) = [-3 1 1] 5 Finite differences scheme Centered Upwind Upwind-Biased Number of points Min. Extr. Max. Extr. Min. Extr. Max. Extr. Min. Extr. Max. Extr. 11 −5 5 - - - - 9 −4 4 - - - - 7 −3 3 - - - - 5 −2 2 - - −3 1 4 - - −3 0 −2 1 3 −1 1 −2 0 - - 2 - - −1 0 - - Table 1: Extreme points in the patterns for constructing matrix eq using the different finite difference methods included in the MatMOL toolbox In the Table 1 we show the extreme points in the patterns that should be used in the construction of the eq matrix for the different finite difference schemes included in the MatMOL toolbox: 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 first derivative of the field. v u ne u 1 X (u(i + 1, j) − u(i, j))2 m(i) = tα + ne (z(i + 1) − z(i))2 j=1 • ’der2’. Monitoring based on the second derivative of the field. v u 2 u(i,j) ne ∂ 2 u(i+1,j) u + ∂ ∂z 1 X 2 ∂z 2 t m(i) = α + ne 2 j=1 • ’derf1’. Monitoring based on the first derivative of the flux. v u ne u 1 X (f l[u(i + 1, j)] − f l[u(i, j)])2 m(i) = tα + ne (z(i + 1) − z(i))2 j=1 This choice is easily implemented in Matlab c , for instance: % Monitoring choice choice = ’der1’; 6 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 adaptation 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 files). In the moving grid algorithm, the state variables and the coordinates of the internal (moving) grid points are taken together to solve the problem. The new vector x constructed from the state variables and grid coordinates is of the form: x = [u11 , u21 , ..., un1 e , z1 , u12 , u22 , ..., un2 e , z2 , ..., u1n , u2n , ..., unne , zn ]T where the subindices indicate the position in the grid and the superindices denote the field considered. In this way, uji refers to the value of the field uj in the grid point zi . The next block in the script Main file.m will assemble state variables and grid coordinates 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); 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 = b(x, t) dt 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 x being the global coordinates as indicated in Main file.m. This is, a vector which includes the state variables and the grid coordinates. 7 In this file the user will define the right hand side terms (b(x, t)). These terms are implemented in matrix udot(:,k) where k indicates the number of the equation. In this file, the user will also choose the finite difference schemes for the computation of the spatial derivatives. It should be noted that vector x contains only the interior (moving) grid points, nevertheless, b(x, t) also makes use of the boundary points z0 , zn+1 , uj0 , ujn+1 with j = 1, .., ne . The structure of b(x, t) ∈ R(ne +1)n×1 is: dune du1 du2 dune du1 du2 dune du11 du21 , , ..., 1 , g1 , 2 , 2 , ..., 2 , g2 , ..., n , n , ..., n , gn b(x, t) = dt dt dt dt dt dt dt dt dt T where µ 1 + 2µ µ 1 µ 1 + 2µ µ 1 − + − − − + − gi = Mi ∆zi+1 ∆zi ∆zi−1 Mi−1 ∆zi ∆zi−1 ∆zi−2 For details see [1, 3]. In the right vect.m function the user will first 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 The second block will make a call to the function which implements the boundary conditions (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 define the nonlinear functions, choose the finite difference scheme and compute the spatial derivatives of the fields and nonlinear functions: % Nonlinear functions f1 = u(:,3); f2 = (gam-1)*u(:,3)- ((gam-3)*u(:,2).ˆ2)./(2*u(:,3)); 8 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) % 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 specified in Dirichlet form. This means that for Neumann and Robin type a transformation have to be performed. The first 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: 9 • To compute the value of the field in the boundary by hand. If, for instance, we have the following boundary condition at the right boundary: a1 ∂uk (zR ) + a0 uk (zR ) = f (t) ∂z (8) we can approximate the first derivative by a simple finite difference scheme so that ∂uk (zR ) u(n + 2, k) − u(n + 1, k) = ∂z z(n + 2) − z(n + 1) Substituting this expression into equation (8) and after some algebraic manipulations we obtain: u(n + 1, k) z(n + 2) − z(n + 1) a1 a0 + z(n + 2) − z(n + 1) f (t) + a1 u(n + 2, k) = • 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. % 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 10 % 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 finite 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 convective. 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 specified. It is necessary, then, to intoduce a “fictitious” boundary condition. From the point of view of time integration, the most flexible way of defining such conditions is to make the first spatial derivative of the corresponding function to be constant. If we consider the first point of the spatial domain, this translates into: ∂u(2, k) ∂u(1, k) = ∂z ∂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
Source Exif Data:
File Type : PDF File Type Extension : pdf MIME Type : application/pdf PDF Version : 1.4 Linearized : No Page Count : 11 Producer : pdfTeX-1.40.3 Creator : TeX Create Date : 2009:10:21 08:43:11+02:00 Modify Date : 2009:10:21 08:43:11+02:00 Trapped : False PTEX Fullbanner : This is pdfTeX using libpoppler, Version 3.141592-1.40.3-2.2 (Web2C 7.5.6) kpathsea version 3.5.6EXIF Metadata provided by EXIF.tools