http://www.springer.com/354022842X

2

Implementationofmodelproblems

In this chapter we describe the implementation of two stationary model prob-lems (the linear Poisson equation and a nonlinear reaction–diﬀusion equation) and of one time dependent model problem (the heat equation). Here we give an overview how to set up anALBERTAprogram for various applications. We do not go into detail when refering toALBERTAdata structures and functions. A detailed description can be found in Chapter 3. We start with the easy and straight forward implementation of the Poisson problem to learn about the basics ofALBERTA. The examples with the implementation of the nonlinear reaction-diﬀusion problem and the time dependent heat equation are more involved and show the tools ofALBERTAfor attacking more complex prob-A lems. Removing all L T X descriptions of functions and variables results in the E source code for the adaptive solvers. During the installation ofALBERTA(described in Section 2.4) a subdi-rectoryDEMOwith sources and makeﬁles for these model problems is cre-ated. The corresponding ready-to-run programs can be found in the ﬁles ellipt.c,heat.c, andnonlin.c,nlprob.c,nlsolve.cin the subdirec-toryDEMO/src/Common/. Executable programs for diﬀerent space dimensions can be generated in the subdirectoriesDEMO/src/1d/,DEMO/src/2d/, and DEMO/src/3d/by callingmake ellipt,make nonlin, andmake heat. Graph-ics output for all problems is generated via a routine void graphics(MESH *mesh, DOF_REAL_VEC *u_h, REAL (*get_est)(EL *)); which shows the geometry given bymesh, as well as ﬁnite element func-tion values given byu h, or local estimator values when parameterget est is given, all in separate graphics windows. The source of this routine is DEMO/src/Common/graphics.c, which is self explaining and not described here in detail.

56

2 Implementation of model problems

2.1 Poisson equation

In this section we describe a model implementation for the Poisson equation d −Δu=finΩ⊂R, u=gon∂Ω.

This is the most simple elliptic problem (yet very important in many appli-cations), but the program presents all major ingredients for general scalar stationary problems. Modiﬁcations needed for a nonlinear problem are pre-sented in Section 2.2.

Fig. 2.1.Solution of the linear Poisson problem and corresponding mesh. The pictures were produced by GRAPE.

Data and parameters described below lead in 2d to the solution and mesh shown in Fig. 2.1. The implementation of the Poisson problem is split into several major steps which are now described in detail.

2.1.1 Include ﬁle and global variables

AllALBERTAsource ﬁles must include the header ﬁlealberta.hwith all ALBERTAtype deﬁnitions, function prototypes and macro deﬁnitions: #include <alberta.h> For the linear scalar elliptic problem we use four global pointers to data struc-tures holding the ﬁnite element space and components of the linear system of equations. These are used in diﬀerent subroutines where such information cannot be passed via parameters. static const FE_SPACE *fe_space; static DOF_REAL_VEC *u_h = nil; static DOF_REAL_VEC *f_h = nil; static DOF_MATRIX *matrix = nil;

2.1 Poisson equation

57

fe spacepointer to the actually used ﬁnite element space; it is initialized: a by the functioninit dof admin()which is called byGET MESH(), see Section 2.1.4; u h: a pointer to a DOF vector storing the coeﬃcients of the discrete so-lution; it is initialized on the ﬁrst call ofbuild()which is called by adapt method stat(), see Section 2.1.7; f h: a pointer to a DOF vector storing the load vector; it is initialized on the ﬁrst call ofbuild(); matrix: a pointer to a DOF matrix storing the system matrix; it is initialized on the ﬁrst call ofbuild(); The data structureFE SPACEis explained in Section 3.2.14,DOF REAL VECin Section 3.3.2, andDOF MATRIXin Section 3.3.4. Details about DOF adminis-trationDOF ADMINcan be found in Section 3.3.1 and about the data structure MESHfor a ﬁnite element mesh in Section 3.6.1.

2.1.2 The main program for the Poisson equation

The main program is very simple, it just includes the main steps needed to implement any stationary problem. Special problem-dependent aspects are hidden in other subroutines described below. We ﬁrst read a parameter ﬁle (indicating which data, algorithms, and solvers should be used; the ﬁle is described below in Section 2.1.3). Then we initialize the mesh (the basic geometric data structure), and read the macro triangulation (including an initial global reﬁnement, if necessary). The subdirectoriesMACROin theDEMO/src/*ddirectories contain data for several sample macro triangulations. How to read and write macro triangula-tion ﬁles is explained in Section 3.2.16. The macro ﬁle name and the number of global reﬁnements are given in the parameter ﬁle. Now, the domain’s ge-ometry is deﬁned, and a ﬁnite element space is automatically generated via theinit dof admin()routine described in Section 2.1.4 below. A call to graphics()displays the initial mesh. The basic algorithmic data structureADAPT STATintroduced in Section 3.13.1 speciﬁes the behaviour of the adaptive ﬁnite element method for sta-tionary problems. A pre–initialized data structure is accessed by the func-tionget adapt stat(); the most important members (adapt>tolerance, adapt>strategy, etc.) are automatically initialized with values from the pa-rameter ﬁle; other members can be also initialized by adding similar lines for these members to the parameter ﬁle (compare Section 3.13.4). Eventu-ally, function pointers for the problem dependent routines have to be set (estimate,get el est,build,solve). Since the assemblage is done in one step after all mesh modiﬁcations, onlyadapt>build after coarsenis used, no assemblage is done before reﬁnement or before coarsening. These addi-tional assemblage steps are possible and may be needed in a more general application, for details see Section 3.13.1.

58

2 Implementation of model problems

The adaptive procedure is started by a call ofadapt method stat(). This automatically solves the discrete problem, computes the error estimate, and reﬁnes the mesh until the given tolerance is met, or the maximal number of iterations is reached, compare Section 3.13.1. Finally,WAIT REALLYallows an inspection of the ﬁnal solution by preventing a direct program exit with closure of the graphics windows.

int main(int argc, char **argv) { FUNCNAME("main"); MESH *mesh; int n_refine = 0; static ADAPT_STAT *adapt; char filename[100];

/**/ /* first of all, init parameters of the init file */ /**/

init_parameters(0,

"INIT/ellipt.dat");

/**/ /* get a mesh, and read the macro triangulation from file */ /**/

mesh = GET_MESH("ALBERTA mesh", init_dof_admin, init_leaf_data); GET_PARAMETER(1, "macro file name", "%s", filename); read_macro(mesh, filename, nil); GET_PARAMETER(1, "global refinements", "%d", &n_refine); global_refine(mesh, n_refine*DIM); graphics(mesh, nil, nil);

/**/ /* init adapt structure and start adaptive method */ /**/

adapt = get_adapt_stat("ellipt", "adapt", 2, nil); adapt>estimate = estimate; adapt>get_el_est = get_el_est; adapt>build_after_coarsen = build; adapt>solve = solve;

adapt_method_stat(mesh, adapt); WAIT_REALLY; return(0); }

2.1 Poisson equation

2.1.3 The parameter ﬁle for the Poisson equation

59

The following parameter ﬁleINIT/ellipt.datis used for theellipt.cpro-gram: macro file name: Macro/macro.amc global refinements: 0 polynomial degree: 3

% graphic windows: solution, estimate, and mesh if size > 0 graphic windows: 300 300 300 % for graphics you can specify the range for the values of % discrete solution for displaying: min max % automatical scaling by display routine if min >= max graphic range: 0.0 1.0

solver: 2 % 1: BICGSTAB 2: CG 3: GMRES 4: ODIR 5: ORES solver max iteration: 1000 solver restart: 10 % only used for GMRES solver tolerance: 1.e8 solver info: 2 solver precon: 2 % 0: no precon 1: diag precon % 2: HB precon 3: BPX precon

error norm: estimator C0: estimator C1: estimator C2:

adapt>strategy: adapt>tolerance: adapt>MS_gamma: adapt>max_iteration: adapt>info:

1 % 1: H1_NORM, 2: L2_NORM 0.1 % constant of element residual 0.1 % constant of jump residual 0.0 % constant of coarsening estimate

2 % 0: no adaption 1: GR 2: MS 3: ES 4:GERS 1.e4 0.5 20 8

WAIT: 1 The ﬁleMacro/macro.amcstoring data about the macro triangulation forΩ= d (0,1) can be found in Section 3.2.16 for 2d and 3d. Thepolynomial degree parameter selects the third order ﬁnite elements. Bygraphic windows, the number and sizes of graphics output windows are selected. This line is used by thegraphics()routine. For 1d and 2d graphics, the range of function values might be speciﬁed (used for graph coloring and height). The solver for the linear system of equations is selected (here: the conjugate gradient solver), and corresponding parameters like preconditioner and tolerance. Parameters for the error estimator include values of diﬀerent constants and selection of 1 2 the error norm to be estimated (H- orL-norm, selection leads to multipli-cation with diﬀerent powers of the local mesh size in the error indicators), see Section 3.14.1. An error tolerance and selection of a marking strategy with

60

2 Implementation of model problems

corresponding parameters are main data given to the adaptive method. Fi-nally, theWAITparameter speciﬁes whether the program should wait for user interaction at additional breakpoints, whenever aWAITstatement is executed as in the routinegraphics()for instance. The solution and corresponding mesh in 2d for the above parameters are shown in Fig. 2.1. As optimal parameter sets might diﬀer for diﬀerent space dimensions, separate parameter ﬁles exist in1d/INIT/,2d/INIT/, and 3d/INIT/.

2.1.4 Initialization of the ﬁnite element space

During the initialization of the mesh byGET MESH()in the main program, we have to specify all DOFs which we want to use during the simulation on the mesh. The initialization of the DOFs is implemented in the function init dof admin()which is called byGET MESH(). For details we refer to Sec-tions 3.2.15 and 3.6.2. For the scalar elliptic problem we need one ﬁnite element space for the discretization. In this example, we use Lagrange elements and we initialize the degree of the elements via a parameter. The correspondingfe spaceis accessed byget fe space()which automatically stores at the mesh informa-tion about the DOFs used by this ﬁnite element space. It is possible to access several ﬁnite element spaces inside this function, for instance in a mixed ﬁnite element method, compare Section 3.6.2.

void init_dof_admin(MESH *mesh) { FUNCNAME("init_dof_admin"); int degree = 1; const BAS_FCTS *lagrange;

GET_PARAMETER(1, "polynomial degree", "%d", °ree); lagrange = get_lagrange(degree); TEST_EXIT(lagrange)("no lagrange BAS_FCTS\n"); fe_space = get_fe_space(mesh, lagrange>name, nil, lagrange); return; }

2.1.5 Functions for leaf data

As explained in Section 3.2.12, we can “hide” information which is only needed on a leaf element at the pointer to the second child. Such information, which we use here, is the local error indicator on an element. For this elliptic problem we need oneREALfor storing this element indicator. During mesh initialization byGET MESH()in the main program, we have to give information about the size of leaf data to be stored and how to transform

2.1 Poisson equation

61

leaf data from parent to children during reﬁnement and vice versa during coarsening. The functioninit leaf data()initializes the leaf data used for this problem and is called byGET MESH(). Here, leaf data is one structure struct ellipt leaf dataand no transformation during mesh modiﬁcations is needed. The details of theLEAF DATA INFOdata structure are stated in Section 3.2.12. The error estimation is done by the library functionellipt est(), see Section 3.14.1. Forellipt est(), we need a function which gives read and write access to the local element error, and for the marking function of the adaptive procedure, we need a function which returns the local error indicator, see Section 3.13.1. The indicator is stored as theREALmemberestimate ofleaf datastruct ellipt and the functionrw el est()returns for each element a pointer to this member. The functionget el est()returns the value stored at that member for each element.

struct ellipt_leaf_data { REAL estimate; };

/*

one real for the estimate

*/

void init_leaf_data(LEAF_DATA_INFO *leaf_data_info) { leaf_data_info>leaf_data_size = sizeof(struct ellipt_leaf_data); leaf_data_info>coarsen_leaf_data = nil; /* no transformation */ leaf_data_info>refine_leaf_data = nil; /* no transformation */ return; }

static REAL *rw_el_est(EL *el) { if (IS_LEAF_EL(el)) return(&((struct ellipt_leaf_data else return(nil); }

static REAL get_el_est(EL *el) { if (IS_LEAF_EL(el)) return(((struct ellipt_leaf_data else return(0.0); }

*)LEAF_DATA(el))>estimate);

*)LEAF_DATA(el))>estimate);

62

2 Implementation of model problems

2.1.6 Data of the diﬀerential equation

Data for the Poisson problem are the right hand sidefand boundary values g. For test purposes it is convenient to have access to an exact solution of the problem. In this example we use the function

as exact solution, resulting in

2 −10|x| u(x) = e

2 −10|x| ∇u(x) =−20xe

and 2 2−10|x| f(x) =−Δu(x) =−(400|x| −20d) e. d Here,ddenotes the space dimension,Ω⊂R. The functionsuandgrd uare the implementation ofuand∇uand are optional (and usually not known for a general problem). The functionsgandfare implementations of the boundary values and the right hand side and are not optional.

static REAL u(const REAL_D x) { return(exp(10.0*SCP_DOW(x,x))); }

static const REAL *grd_u(const REAL_D x) { static REAL_D grd; REAL ux = exp(10.0*SCP_DOW(x,x)); int n;

for (n = 0; n < DIM_OF_WORLD; n++) grd[n] = 20.0*x[n]*ux;

return(grd); }

static REAL g(const REAL_D x) { return(u(x)); }

/* boundary values, not optional */

static REAL f(const REAL_D x) /* Delta u, not optional { REAL r2 = SCP_DOW(x,x), ux = exp(10.0*r2); return((400.0*r2 20.0*DIM)*ux); }

*/

2.1 Poisson equation

2.1.7 The assemblage of the discrete system

63

For the assemblage of the discrete system we use the tools described in Sections 3.12.2, 3.12.4, and 3.12.5. For the matrix assemblage we have to provide an element-wise description of the diﬀerential operator. Following the description in Section 1.4.7 we provide the functioninit element()for an initialization of the operator on an element and the functionLALt()for the computation t of det|DFS|ΛAΛon the actual element, whereΛis the Jacobian of the barycentric coordinates,DFSthe the Jacobian of the element parameteriza-tion, andAthe matrix of the second order term. For−Δ, we haveA=idand t det|DFS|ΛΛis the description of diﬀerential operator since no lower order terms are involved. For passing information about the JacobianΛof the barycentric coordi-nates and det|DFS|from the functioninit element()to the functionLALt() we use the data structureinfostruct op which stores the Jacobian and the determinant. The functioninit element()calculates the Jacobian and the determinant by the library functionel grd lambda()and the functionLALt() t uses these values in order to compute det|DFS|ΛΛ. Pointers to these functions and to one structureinfostruct op are mem-bers of a structureOPERATOR INFOwhich is used for the initialization of a function for the automatic assemblage of the global system matrix. For more general equations with lower order terms, additional functionsLb0,Lb1, and/orchave to be deﬁned at that point. The full description of the func-tionfill matrix info()for general diﬀerential operators is given in Sec-tion 3.12.2. Currently, the functionsinit element()andLALt()only work properly for== DIMDIM OF WORLD . The initialization of theEL MATRIX INFO structure is only done on the ﬁrst call of the functionbuild()which is called byadapt method stat()during the adaptive cycle (compare Section 3.13.1). By callingdof compress(), unused DOF indices are removed such that the valid DOF indices are consecutive in their range. This guarantees op-timal performance of the BLAS1 routines used in the iterative solvers and admin>size usedis the dimension of the current ﬁnite element space. This dimension is printed for information. On the ﬁrst call,build()also allocates the DOF vectorsu handf h, and the DOF matrixmatrix. The vectoru hadditionally is initialized with zeros and the function pointers for an automatic interpolation dur-ing reﬁnement and coarsening are adjusted to the predeﬁned functions in fe space>bas fcts. The load vectorf hand the system matrixmatrixare newly assembled on each call ofbuild(). Thus, there is no need for interpo-lation during mesh modiﬁcations or initialization. On each call ofbuild()the matrix is assembled by ﬁrst clearing the matrix using the functionclear dof matrix()and then adding element con-tributions byupdate matrix(). This function will callinit element()and LALt()on each element.