Home

DROPS-Package User's guide - Institut für Geometrie und

1. 1 0 M tn G t Bl p 4 36 Assume that is such that Bu 0 The sequence ti n 9 defined by the 6 scheme 4 35 satisfies 4 36 and also Bu 0 for all n We use p as a Lagrange multiplier to enforce B 0 as follows Given i with Bu 0 define B S to BM to G to and for n gt 0 let d p be such that url d AP OM int GE tng BTB 1 Min Clin BTP Bu t 0 4 37 holds Note that in this saddle point type system the projection Q is not used Due to 4 36 this system has a solution If we assume that for each n the saddle point problem 4 37 has a unique solution which is true for At sufficiently small then this yields the solution of the 0 scheme in 4 35 Remark 5 If the mass matrix M does not depend on t then we can take p s BM 0G tg 1 6 G a t 1 k 1 and instead of 4 36 we obtain grt d At which then results in the standard 6 scheme for a one phase Navier Stokes equation as described in section 3 3 1 M7 6G 1 1 1 1 0 G d tr M BT pt For 6 0 the method in 4 37 corresponds to the explicit Euler method applied to 4 34 which is not very useful do to its poor stability properties We consider 0 Z 0 In this case in 4 37 inverses of both M t and M t occur The latter can be avoided by introducing an additional variable leading to a more convenient but equivalent formulation of
2. 7 5 5 ii f E c EM 2 E K 5 Ji 3 26 With A zz Q4 S w S w V w and rf f Axt BT y we obtain the nonlinear iterative method x oxky Qi Q4 B wv B Qji xt f2 E A 3 27 y E y V B Qjri x f2 Thus we obtain an inexact Uzawa method with the following algorithmic structure x 0 A VO T0 y a given starting vector rji fi Ax B y fork 0 wie x Q rt 3 28 z V Bw fi sepe ee Q B z yH yz riti rf AGT x pig Here W is as follows Result of 2 PCG iterations with starting vector 0 ace oe preconditioner Qs applied to Sv Bw fo 3 29 This algorithm consists of an inner outer iteration An analysis of this method is given in 10 where it is shown that in the inner iteration 3 29 one should use very small values 1 2 and that this inexact Uzawa method is an efficient solver for the saddle point problem 3 20 provided we have good preconditioners Q4 of A and Qs of S 3 30 For the preconditioner Q4 we use a standard symmetric multigrid V cycle Gauss Seidel smoother An appropriate preconditioner Qs is due to Cahouet and Chabard 36 We briefly explain this method For further information and a convergence analysis we refer to the literature 10 For g L Q consider the Neumann problem find w H Q n Q such that Vw V g forall de H Q nQq 3 31 Let T be the stiffness matrix of the Galerkin discretization of this
3. 9 1 DATA DISTRIBUTION 81 Figure 9 2 Domain decomposition for P 8 processors Figure 9 2 shows a domain decomposition for P 8 processors A domain decomposition of level 7 automatically induces a distribution of the numerical data on level j Without loss of generality we assume that the global numbering 7 1 N is associated with the finest level J Let J p 1 Np be a local numbering of the degrees of freedom of the local triangulation 77 p on processor p 1 p P Then the relation between a local number i 7 p and its global counterpart j J is given by the coincidence matrix Ip E RN 1 if degree of freedom with global number j 7 exists Lie on processor p with local number i 7 p 0 else Degrees of freedom which are located on multiple processors form the so called processor bound ary Definition 9 Accumulated and distributed storage For a global vector Z R the sequence Za IZ Ipi e RM x x RAP is called the corresponding accumulated vector That means that for unknowns on a processor boundary each adjacent processor stores the same global value The sequence q 1 p of vectors Zp IR is called the distributed vector corre sponding to 7 if ge I p 1 In this case the global value of an unknown on a processor boundary is the sum of all local values stored on the adjacent processors The same holds for entries of a distributed matr
4. The temporal and spatial discretization of the level set equation are not mass conserving Due to the surface tension we usually lose mass from Q1 This loss of mass is reduced if the grid is refined Such finer grids however result in higher computational costs Therefore we introduce another strategy to compensate for the mass loss After each time step we shift the interface in normal direction such that the volume of Q at current time is the same as at time t 0 To realize this we exploit the fact that the level set function is close to a signed distance function In order to shift the interface over a distance d in outward normal direction we only have to subtract d from the level set function Let V 9 vol x d x lt 0 denote the volume of corresponding to a level set function and let be the discrete level set function at a given time We have to find d c R such that V n d vol Q9 0 0 holds In order to keep the number of evaluations of V low we use a method with a high rate of convergence namely the Anderson Bj rk method 17 to solve this equation We then set Gy bn d and discard yp 4 3 TIME INTEGRATION 43 Note that this strategy only works if Q consists of a single component If there are multiple components mass must be preserved for each of them In this case the algorithm can be modified to shift only locally Discontinuities that may occur in the level set functi
5. These parameters are obtained easily using an object of the class ParamMesszelleCL file Level set params h and can be accessed through its members Note that the names of the param eters in the parameter file and the corresponding class members are not identical 12 2 2 Structure of the program The following libraries are needed include geom multigrid h include geom builder h include navstokes instatnavstokes2phase h include stokes integrTime h include num stokessolver h include out output h include out ensightOut h include levelset coupling h include levelset params h include levelset adaptriang h include num bndData h include lt fstream gt Similar to the example in chapter 11 this program consists of e the coefficient class ZeroFlowCL as the template parameter for the problem class Instat NavierStokes2PhaseP2P1CL e the function Null for the initial and Dirichlet boundary conditions of the velocity e the function DistanceFct to initialize the levelset function e the function sigmaf which defines the surface tension If the surface tension is variable one needs an additional function for its gradient e the functions Strategy and main 12 2 IMPLEMENTATION 105 12 2 3 The function main In the function main we read the input parameters from the the parameter file into a global object C of the class ParamMesszelleCL std ifstream param argv 1 param gt gt C para
6. 1 0 q i a 4 44 cy ye F 0 We reformulate this method as follows We decompose d as qd Z w r T We only consider 0 4 0 From C y 0q 9F y 0 and 9gq we get for n gt 0 AEH y 1 8 gt n 1 ur tiet M E gt n 1 gt n E OF 2 t E a 1 E t bw 0 F3 y 6r gH M9 MER and on l in gar I Aj D 0 gt n 1 gt n gwrt LA 1 6 w In this form we have exactly the same scheme as in 4 41 For n 0 we need starting values 4 a gt gt 0 Z d to W to From E u H we get E w H and thus the same starting value for W as in 4 40 For z we consider 4 31 and obtain dB d Be N and thus with the same notation as in the previous section SHB B t M t G d t Fan 48 CHAPTER 4 NAVIER STOKES EQUATIONS FOR TWO PHASE FLOW NS2 From this we obtain the system 5 30 49 72 22 dB S to p B to M 6 G u 6 g fr T Note that if B does not depend on t this system is the same as the one in 4 40 Given this p the starting value Z can be determined from 20 ao 50 M z zu G amp g fr Bip to W to 4 45 which is the same as in 4 40 We see that if B is independent of t then we obtain the same starting values as in 4 40 and hence the method based on the DAE system 4 43 is exactly the same as the one derived in the previous section 4
7. Pw if Av T Tm mini lt j lt m v Er Q otherwise Since a vertex is part of multiple tetrahedra we still have to decide which tetrahedron T we use to determine the approximate distance d v d v min dr v v e V T and T Tr 4 27 All vertices v V T T Tr now have values d v Far field phase In this second phase of the FMM we extend the approximate distance function d v to the neighbors of those vertices which are already assigned a value To organize this procedure we introduce three sets FAR to contain those vertices which are not affected by the algorithm yet FIN containing those that are already assigned a final value and ACT which contains the vertices that are currently updated i e those adjacent to FIN Let V Q denote the set of all vertices in and 7 v denote the set of tetrahedra adjacent to v The three sets are initialized as follows FIN veV 9 v is assigned a value d v during initialization ACT veV 9 v has a neighbor in FIN FAR v V Q v has no neighbor in FIN We determine values d v for all v ACT in a similar way as in the initialization phase Take v ACT T T v such that V T n FIN There are three possible cases namely IV T n FIN 1 2 or 3 If V T n FIN 1 say V T n FIN w we define dr v d w v w 42 CHAPTER 4 NAVIER STOKES EQUATIONS FOR TWO PHASE FLOW NS2 If T has m 2 or m 3 vertices co
8. ini c 1 2 C cD D O C diffusityF D 1 C diffusityD H C equilibrium jump e In the function Strategy we define an object c of the TransportP1CL class with cBndDataCL Bnd c 6 c bc c bfun TransportP1CL c MG Bnd c Bnd v theta 1 D H amp Stokes v lset t 0 C dt C transp iter C transp tol IdxDescCL cidx amp c idx c CreateNumbering MG GetLastLevel cidx c ct SetIdx cidx c Init amp Initialcneg amp Initialcpos c Update c SetTimeStep C dt The solutions of the transport problem can be exported to the ensight output files const string datc filename c datct filenamet ct ensight DescribeScalar Concentration datc true ensight DescribeScalar TransConc datct true ensight putScalar datc c GetSolution 0 ensight putScalar datct c GetSolution c ct 0 In each time step after solving the two phase flow part we solve the transport problem with the function DoStep of c If the grid is modified the concentration should be updated with the member functions UpdateTriang of adap using the address of c as the third parameter and Update of c for int step 1 step lt C num_steps step 1 Std cerr lt lt Schritt lt lt step lt lt t t lt lt step C dt lt lt n cpl DoStep C cpl_iter 13 3 RESULTS 113 c DoStep step C dt lset GetInfo maxGradPhi Volume bary drop min drop max drop infofile lt lt Stokes
9. supp R V supp v lt 1 It is clear that if functions with very small support are deleted from the space WI this will not influence the approximation quality of the XFEM space WT significantly Therefore we introduce a smaller space in which basis functions from WF with very small support are deleted Avoiding very small supports has advantages for example if the contributions are dominated by rounding errors We will explain how we chose the maximal size of these small supports in order to maintain optimal approximation properties of the resulting reduced XFEM space Let a gt 0 gt 0 be given parameters Let Z7 C Jr be the index set such that for all k Jr RE ne lt h forall T C supp ux n Nr 4 23 kIT Remark 4 Note that for a function Riv WI k TF we have R vyllimr lWellarne for i all T 7 Furthermore because v li r ch 5 for 0 1 the condition 4 23 can be replaced by the following one Wella uh for all T C supp vx n Qr 4 24 The constant may differ from in 4 23 We define the reduced spaces W C WF by W span Rid k ETT i 12 and a modified XFEM space Wl W W amp W7 For this space the following approximation property holds cf 13 42 SPATIAL DISCRETIZATION 39 Theorem 3 We assume 7 1450 to be quasi uniform For 0 l m 2 the following holds m Iu vll uR e t ip ula Quos for all u H N UM vc From t
10. 0 We list a few important properties of this method e Both the input and output are multilevel triangulations The method yields stable and consistent triangulations Local refinement and coarsening are treated in a similar way The implementation uses only the hierarchical decomposition of M This allows relatively simple data structures without storage over head The costs are proportional to the number of tetrahedra in 77 For a detailed discussion of these and other properties we refer to the literature 22 27 63 1 In our implementation we use the multilevel refinement algorithm described in 1 Chapter 3 Navier Stokes equations for one phase flow NS1 3 1 Weak formulation We consider the Navier Stokes model as in 1 1 For simplicity we take p 1 u constant and homogeneous Dirichlet boundary conditions for u on Q For the weak formulation of this problem we introduce the spaces vH Q RO teer a 0 3 1 A weak formulation of this problem is as follows Determine u t V and p t Q such that m u t v a u t v c u t u t v b v p t b u t a for almost all t 0 T Here and in the remainder denotes the component wise L scalar product In this variational formulation we used the bilinear forms lan vac se Es on A Vu Vv de b u g f adivuar Q and the trilinear form c u v w nc Vv w dz Q 3 2 Spatial discretization 3 2 1 Galerki
11. 4 37 This is done as follows Define Z M t G u t B p k2 0 4 3 TIME INTEGRATION 45 Le M t9 Z G u9 to B p og u _ ur 1 9 2 gt 0 A 1 8 z k gt 0 Using thus 4 37 can be reformulated as grt _ qn M ta r 0B p OGG tay 1 AM tn 2 Bu t 0 4 38 on l 3n garth 8z for n gt 0 and a starting value Z as defined above Application of the 0 scheme to the level set equation 4 9 results in 2n F1 2n Sy gt n antl zr gd 0B a1 Bat 1 oE HRG This can be reformulated using a new variable w E f H u which satisfies for 0 4 0 2n 1 on wir 1 0 w m At eee resulting in gt n 1 g 1 Eat IH ET 1 6 e a twn 4 39 Combining these results and inserting the notation for G we obtain for 0 Z 0 the following Sia gt 5 zn 1 coupled nonlinear system for p gt p i d 46 CHAPTER 4 NAVIER STOKES EQUATIONS FOR TWO PHASE FLOW NS2 Given uf di determine Z and w as follows mB E w w For n gt 0 n 1 2ndl u s 2n41 gt 3n4l M T ILN amp u Du n l g fr 4B pr 5n 4 1 1 0 M Z gt n 1 gt n Ba 0H u 9 2 gt n 1 n 1l on E u IR 1 60 E u w on l d At H 6 1 8 w n 1 0 w Remark 6 The derivation above shows that the scheme in 4 40 4 4
12. 5 4 The diffusion coefficient is piecewise constant D D Da D1 H In the interface condition we use the notation c for cjg restricted to the interface The constant Cg gt 0 is given Henry s constant The model has to be combined with suitable initial and boundary conditions Note that the transport equation needs as input the flow field u that results from the Navier Stokes equations If the surface tension coefficient 7 is assumed to be independent of the concen tration c there is no dependence of the Navier Stokes equation 5 1 on the transport problem 5 2 Due to this we can use a simple one direction decoupling between the Navier Stokes equation and the transport problem as described in section 5 4 The weak formulation of the Navier Stokes equations for two phase flow is given in section 4 1 Below we describe a weak formulation of the convection diffusion equation 5 2 Note that due to the Henry interface condition in 5 4 the concentration c is discontinuous across the inter face for Cy Z 1 To eliminate the discontinuity we introduce a piecewise constant weighting 51 52CHAPTER 5 TWO PHASE FLOW WITH TRANSPORT OF A DISSOLVED SPECIES NS24 M function 1 in Qi i e if z t 0 Cn 9 t f Cy in ie if d z t gt 0 We define the scaled function Cy c In view of 5 4 this function should be continuous maybe in some weak sense across the interface For this transforme
13. For example the iterative linearization by a fixed point method of the one phase Navier Stokes equations as described in section 3 4 results in a linear saddle point problem Oseen equation For this linear Oseen problem one can use a preconditioned MINRES method In the preconditioner an inner iterative solver for example multigrid may be used Thus in this example one has three nested iterative methods a fixed point linearization method an Oseen solver and an inner preconditioner In view of this we use a template mechanism to specify the inner solution components as template parameters This enables an easy plug in of different solution components to test and compare reasonable combinations of solvers available from the DROPS solver toolbox Furthermore this technique assures efficient code since the compiler can perform full code optimization for the template specialization which is known at the moment of compilation Example 8 1 As an illustrative example for the template plug in mechanism we give a piece of code for the definition of a Stokes solver cf section 3 5 2 preconditioner for upper left block preconditioner typedef SSORPcCL ULPcPcT ULPcPcT ULPcPc preconditioner for upper left block typedef PCGSolverCL lt ULPcPcT gt ULSolverT ULSolverT ULsolver ULPcPc typedef SolverAsPreCL lt ULSolverT gt ULPcT ULPcT ULPc ULsolver Schur complement preconditioner typedef ISPreCL SchurPcT SchurPcT Sch
14. H u n Culo u Vj dz A EREK Ad I D V Vii dz Determine c t such that c 0 is given and Mo H for all t 0 T To determine the entries of the matrices M H and A one has to compute integrals over tetrahe dra in which there are discontinuities in the coefficient functions Cy and D across the interface which is approximated by T For this the methods discussed in section 4 2 5 are used 5 2 2 Nitsche s method combined with XFEM In this section we use a technique from 55 that results in a finite element discretization method with an optimal order of convergence Forthcoming 5 3 Time integration We recall the 0 scheme from section 3 3 1 urew _ yold At 1 8 F u H OF um 0 0 1 54CHAPTER 5 TWO PHASE FLOW WITH TRANSPORT OF A DISSOLVED SPECIES NS24 M If this method is applied to the discrete transport equation in 5 6 with F E M E A 9 0 H t u 2 6 we obtain the following fully discrete problem with t u u t ert _ e A 6 M G 7 A 9 He u e AtoM 1 A zs H t u art This can be rewritten as M At0 A t 4 H t utes M 1 At 1 BMA A 9 HE u e Thus per time step we have to solve two linear systems if 0 0 1 and only one linear system if 0 1 5 4 Decoupling and linearization As can be seen from the model 5 1 5 4 there is only a coup
15. It contains an array of BndSegDataCL lt BndValT gt objects one for each boundary segment X cf Section 8 1 Each BndSegDataCL object stores the boundary condition and a function pointer for evaluating the corresponding boundary values of type BndValT The choice of the template parameter BndValT depends on whether the boundary condition applies to a scalar double or vector valued Point3DCL quantity The boundary condition of type BndCondT can be one of e DirBC DirOBC for inhomogeneous and homogeneous Dirichlet boundary conditions re spectively e NatBC NatOBC for inhomogeneous and homogeneous natural boundary conditions respec tively 8 5 TOOLS FOR SPATIAL DISCRETIZATION 71 e Per1BC Per2BC for periodic boundary conditions denoting corresponding boundaries WallBC and OutflowBC are alias names for DirOBC and NatOBC respectively Coefficients As an example of how the coefficients of a specific partial differential equation sre represented in the implementation we consider a scalar convection diffusion problem for the unknown function u u x t ut v z t Vu div a x t Vu f x t in Q x to ty This type of problem is defined by the problem class InstatPoissonP1CL The corresponding PoissonCoeffCL contains the functions v z t a z t and f x t as static member functions For the two phase flow problem 4 1 4 2 the corresponding coefficient class stores quantities such as densities p and dynamic viscosit
16. Phenomena with Moving Boundaries pp 134 146 VDI Verlag D sseldorf 2004 BIBLIOGRAPHY 119 31 32 33 34 35 36 37 38 46 47 D BOTHE J PRUss G SIMONETT Well posedness of a two phase flow with soluble surfactant In M CHIPOT J ESCHER ed Nonlinear Elliptic and Parabolic Problems pp 37 61 Birkh user 2005 J U BRACKBILL D B KOTHE C ZEMACH A continuum method for modeling surface tension J Comput Phys 100 pp 335 354 1992 D BRAESS Finite Elements Theory Fast Solvers and Applications in Solid Mechanics Cambridge University Press Cambridge 2007 J H BRAMBLE J E PASCIAK Iterative techniques for time dependent Stokes problems Comput Math Appl 33 No 1 2 pp 13 30 1997 M O BRISTEAU R GLOWINSKI J PERIAUX Numerical methods for the Navier Stokes equations Application to the simulation of compressible and incompressible flows Computer Physics Report 6 pp 73 188 1987 Cahouet J Chabard J P Some fast 3D finite element solvers for the generalized Stokes problem Int J Numer Methods Fluids 8 1988 869 895 B CHAMBERLAIN Graph Partitioning Algorithms for Distributing Workloads of Parallel Computations TR 98 10 03 1998 Y C CHANG T Y Hou B MERRIMAN S OSHER A level set formulation of Eulerian interface capturing methods for incompressible fluid flows J Comput Phys 124 pp 449 464 1996 K DECKELNICK G DZIUK Me
17. e a one phase Navier Stokes problem subdirectory NAVSTOKES e a two phase Navier Stokes problem subdirectory LEVELSET As an example go to the subdirectory navstokes cd navstokes Now run the make utility at the shell prompt in the current directory make lt name gt Here lt name gt may be one of the following e nsdrops for an stationary Navier Stokes problem e insdrops for an instationary Navier Stokes problem e insadrops for an instationary Navier Stokes problem with adaptive grid refinement e drivcav for the driven cavity problem described in Chapter 11 For example we continue with the option drivcav make drivcav 10 4 Running DROPS If the make finished successfully the executable file drivcav has been constructed The simplest way to run the program is a command call such as drivcav 1e 7 200 0 1 100 1 200 0 100 0 1 0 1 A detailed description of this example is given in Section 11 2 4 Another way to run a program in DROPS is the usage of the parameter file name param For example to run the program risingBubbleAdap from the drops levelset subdirectory in which the simulation of a rising bubble is implemented we call risingBubbleAdap risingBubbleAdap param A detailed description of this example is given in Chapter 12 96 CHAPTER 10 HOW TO GET STARTED Chapter 11 Navier Stokes equations for a one phase flow 11 1 Introduction We consider a driven cavity problem of the form 98 iE ea
18. pressible flows with surface tension J Comp Phys 224 pp 40 58 2007 M LARIN A REUSKEN A comparative study of efficient iterative solvers for generalized Stokes equations Numer Lin Alg Appl 15 pp 13 34 2008 M A OLSHANSKII A REUSKEN Analysis of a Stokes interface problem Nu mer Math 103 pp 129 149 2006 M A OLSHANSKI J PETERS A REUSKEN Uniform preconditioners for a parameter dependent saddle point problem with application to generalized Stokes interface equations Numer Math 105 No 252 pp 159 191 2006 M A OLSHANSKII A REUSKEN A Stokes interface problem stability finite element analysis and a robust solver In P Neittaanm ki et al Eds Proceedings of European Congress on Computational Methods in Applied Sciences and Engineering ECCOMAS 2004 M A OLSHANSKII A REUSKEN J GRANDE An Eulerian Finite Element method for elliptic equations on moving surfaces IGPM Preprint 2008 Submitted to SIAM J Numer Anal J PETERS V REICHELT A REUSKEN Fast iterative solvers for discrete Stokes equations SIAM J Sci Comput 27 No 2 pp 646 666 2005 A REUSKEN Multilevel techniques for two phase incompressible flows Proceedings of the 8th European Multigrid Conference EMG 2005 2005 A REUSKEN AND V REICHELT Mit Gittern auf der Spur von Tropfen TerraTech Zeitschrift ftir Altlasten und Bodenschutz 6 pp 35 37 2001 A REUSKEN Analysis of an extended pressure finite
19. typedef AdaptFixedPtDefectCorrCL NavStokesCL OseenSolverT gt NSSolverT NSSolverT statsolver NS oseensolver fp_maxiter fp tol deco red The 0 scheme instatsolver is applied for the time integration After that we obtain the reference to the upper left matrix A to set the pointer A of BBTpc InstatNavStokesThetaSchemeCL NavStokesCL NSSolverT instatsolver NS statsolver theta BBTpc SetMatrixA amp instatsolver GetUpperLeftBlock Now we solve the instationary Navier Stokes equations from time begin to time end with the time step timestep In each time step the new time is set via SetTime and the routine DoStep is called for int timestep 1 timestep num timestep timestep 1 double t time begin timestep dt NS SetTime t instatsolver DoStep v p gt Data F After compiling we run the executable with the command of the form drivcav lt fp_tol gt lt fp_maxiter gt lt deco_red gt lt oseen_maxiter gt lt theta gt lt num_timestep gt lt time_begin gt lt time_end gt lt shell_width gt lt c_level gt lt f_level gt In our example we take the following parameters drivcav 1e 7 200 0 1 100 1 200 0 100 0 20 1 11 3 Results The steady state of the flow which is determined by taking a tolerance 1077 in the fixed point iteration is illustrated in Figures 11 1 and 11 2 In these figures the counter rotating eddies in the lower right and left corner can be seen 102 CHAPTER 11 NAVIER ST
20. 335 2001 J GLIMM J W GROVE X L Li K M SHYUE Y ZENG Q ZHANG Three dimensional front tracking SIAM J Sci Comp 19 pp 703 727 1998 J GLIMM M J GRAHAM J GROVE X L Lr T M SMITH D TAN F TANGERMAN Q ZHANG Front tracking in two and three dimensions Comput Math Appl 35 pp 1 11 1998 Greenbaum A Iterative Methods for Solving Linear Systems Frontiers in Applied Math ematics Vol 17 SIAM Philadelphia 1997 S Gross Parallelisierung eines adaptiven Verfahrens zur numerischen L sung partieller Differentialgleichungen Diploma thesis in german 2002 W HACKBUSCH Multi grid Methods and Applications Springer Berlin Heidelberg 1985 W HACKBUSCH Iterative Solution of Large Sparse Systems of Equations Springer Berlin Heidelberg 1994 A HANSBO P HANSBO An unfitted finite element method based on Nitsche s method for elliptic interface problems Comput Methods Appl Mech Engrg 191 pp 5537 5552 2002 C W Hirt B D NICHOLS Volume of Fluid VOF Method for the Dynamics of Free Boundaries J Comp Phys 39 pp 201 225 1981 S HvsiNG A new implicit surface tension implementation for interfacial flows Int J Numer Meth Fluids 2005 A J JAMES J LOWENGRUB A surfactant conserving volume of fluid method for interfa cial flows with insoluble surfactant J Comp Phys 201 2 pp 685 722 2004 GEORGE KARYPIS VIPIN KUMAR A Parallel Algorithm for Multilev
21. 40 4 41 In that case the consistency and stability properties can be derived from the fact that the method is a reformulation of a 6 scheme applied to the system of ODEs 4 34 cf Remark 6 If B depends on then the scheme 4 41 4 40 with a modified starting value for p as in 4 45 results from a formal application of the 0 scheme to the DAE system 4 43 and an analysis of the accuracy and stability properties of this method is not available yet 4 3 8 An implicit Euler type of method We present a variant of an implicit Euler method which is particularly attractive due to its simplicity Starting point is the ODE system for the Navier Stokes part in 4 34 du E E z QOMO Gd t For the discretization the projection operator is treated explicitly whereas the operator G u t which causes the stiffness in this problem is treated semi implicitly If for the latter we apply an implicit Euler method we obtain the discretization grt d At By introducing a Lagrange multiplier p t S t BM t G u t the projection Q t can be eliminated and we obtain the saddle point form Q t M t4 IG TH tn grt At M t4 G um Bp y Bu t 0 4 46 The same idea can be applied to the level set equation and thus we get the following time integration scheme for the coupled problem i 5 0 79 Given amp determine for n gt 0 n gt n u gt n gt n M amp aa a
22. B does not depend on t To obtain a system of ODEs we eliminate the algebraic equation Bu t 0 and the Lagrange multiplier p t as follows In the notation we suppress the dependence on and fr and write M t M t G t G d g fr From BH t 0 here we used that B does not depend on t and substitution of Et from the first equation we obtain S t B t BM t G d t S t BM t BT 4 33 The matrix BT has full rank and thus S t is invertible Using 4 33 we can eliminate p t from the first equation in 4 32 resulting in du 1mTq 1 IQR a I M t B S t B M t G d t Q tM t TIG T t 4 34 44 CHAPTER 4 NAVIER STOKES EQUATIONS FOR TWO PHASE FLOW NS2 The projection Q t J M t BTS t B satisfies BQ t 0 Hence if B 0 0 then the solution t of the ordinary differential equation 4 34 remains in the subspace Ker B To this system of ODEs the 0 scheme is applied resulting in the discretization u At unti en OQ tny M tn 1G d 1 t 1 0 Q t M t G t 4 35 We assume that for each n this system has a unique solution which is the case for At sufficiently small If B 0 then B 0 for all n gt 1 For implementation it is convenient to eliminate the projection Q by introducing a suitable Lagrange multiplier Define p S t BM t 1G t Then 4 35 takes the form grt d Mtn GG tng BTB
23. Math 56 pp 645 666 1990 P BASTIAN Parallele adaptive Mehrgitterverfahren Teubner Stuttgart 1996 P BASTIAN K BIRKEN K JOHANNSEN S LANG N Neuss H RENTZ REICHERT C WIENERS UG A flexible software toolbox for solving partial differential equations Computing and Visualization in Science 1 pp 27 40 1997 P BASTIAN Load balancing for adaptive multigrid methods SIAM J Sci Comput 19 pp 1303 1321 1998 P BASTIAN K BIRKEN K JOHANNSEN S LANG V REICHENBERGER G WITTUM C WROBEL A parallel software platform for solving problems of partial differential equa tions using unstructured grids and adaptive multigrid methods In High performance com puting in science and engineering E Krause and W J ger eds pp 326 339 Springer Berlin 1999 M BEHR Free surface flow simulations in the presence of inclined walls Comput Methods Appl Mech Engrg 191 pp 5467 5483 2002 J Bey Tetrahedral grid refinement Computing 55 pp 355 378 1995 J BEv Finite Volumen und Mehrgitterverfahren f r elliptische Randwertprobleme Ad vances in Numerical Methods Teubner Stuttgart 1998 J BEY Simplicial grid refinement on Freudenthal s algorithm and the optimal number of congruence classes Numer Math 85 pp 1 29 2000 D BorHE M KoEBE H J WARNECKE VOF simulations of the rise behavior of single airbubbles with oxygen transfer to the ambient liquid In F P SCHINDLER ed Transport
24. This rule has seven nodes and corresponding weights as given in table 4 1 nodes in barycentric coordinates A 6 v15 21 B 9 2 15 21 Ag 6 V 15 21 B 9 2 15 21 Table 4 1 Quadrature rule on triangles of degree 5 4 2 6 Re initialization method for the level set function During the evolution the level set function can become distorted and in general loses its property of being an approximate signed distance function At several places in the numerical routines it plays a role that Iva enun gl Go T diss 0x1 0x3 x3 in particular for x close to the interface Thus there is a need for re initialization of the level set function 4 After evolving over a few time steps it is re initialized based on the following criteria we try to fulfil 4 25 the change in the discrete interface I zero level of 5 should be small we aim at mass conservation The condition 4 25 leads to the so called Eikonal equation V 1 in Q dr on Ten For this kind of equation several numerical solution methods are known In particular for re initialization one often introduces a pseudo time variable 7 and considers an instationary problem of the form OY Or with initial condition 0 x old and Sa a regularized sign function Numerical methods can be applied to determine an aproximate stationary solution of this equation which then is taken as re initialization of os S
25. Vol 1e 9 lset Phi Data dphi if C RepFreq amp amp step C RepFreq 0 lset ReparamFastMarching C RepMethod adap UpdateTriang Stokes lset if adap WasModified cpl Update if C VolCorr double dphi lset AdjustVolume Vol 1e 9 lset Phi Data dphi ensight putGeom datgeo step C dt ensight putScalar datpr Stokes GetPrSolution step C dt ensight putVector datvec Stokes GetVelSolution step C dt ensight putScalar datscl lset GetSolution step C dt ensight Commit 108 CHAPTER 12 NAVIER STOKES EQUATIONS FOR A TWO PHASE FLOW In each time step the member function DoStep of cpl is used to solve the levelset and the Navier Stokes equations After that we compute the amount of volume difference in with the routine AdjustVolume and add the correction term to the levelset function The reparametrization is performed after each C RepFreq steps by the member function ReparamFastMarching of levelset The triangulation is then modified and the solutions are interpolated based on the new re initialized levelset function with the function UpdateTriang of adap If the numberings and indices of the vectors are changed the mass matrix for the pressure should also be updated At the end we export the geometry and the solutions to the ensight files with the member functions putGeom putScalar and putVector of ensight 12 3 Results The domain is partitioned into 4 x 20 x 4 cubes a
26. because then the iterative solvers only have to deal with matrices and vectors and not with the grid Then a matrix vector multiplication does not require a loop over all grid entities which saves substantial computational time For the interpretation of a solution vector u however it is necessary to know which vector entries are associated with a certain vertex V for example Here the concept of indices comes into play Index descriptions and numberings For each finite element type used in a solution strategy there exists an associated index An index J is described by an IdxDescCL object It contains the number of degrees of freedom DoF for each simplex type ny ng np nr and the overall number of unknowns Nz To give an example a P index has ny 1 DoF per vertex and ng np np 0 an index for vector valued P5 FE has ny ng 3 DoFs for each vertex and edge and np nr 0 As a next step we have to create a numbering of all degrees of freedom which belong to the index J where degrees of freedom on Dirichlet boundaries are omitted This is done by a function CreateNumbering which is usually a member function of the problem class cf Section 8 4 By this we also obtain the total number of unknowns N 7 which is equal to the dimension of the vectors associated with 7 Thus at the end CreateNumbering sets the value Nz in the corresponding IdxDesc object The numbering is stored by UnknownHandleCL objects contained in the
27. conditions if for the tetrahedral family we take 7 gt o as explained in section 4 2 2 The main new idea of our approach is that for discretization of the problem 6 2 we use a finite element space induced by the continuous linear finite elements on Tj This is done as follows We define a subdomain that contains Tp Wh Ure ST 6 4 This subdomain in IR is connected and partitioned in tetrahedra that form a subset of Tp We introduce the finite element space Vj um vp C wn Us P forall Te Fn 6 5 This space induces the following space on Ij VE vy H T4 3 vy Va vy vulr 6 6 This space is used for a Galerkin discretization of 6 2 determine up Vr with Sr upds 0 such that Vr uj Vr Vy dsp ni frvnds forall yy VI 6 7 INA In with fj a suitable extension of f such that Sr Jndsn 0 Due the Lax Milgram lemma this problem has a unique solution up In 9 we present a discretization error analysis of this method that shows that under reasonable assumptions we have optimal error bounds Remark 8 As far as we know this method for discretization of a partial differential equation on a surface is new We give some comments related to this approach e The family 751550 is shape regular but the family I 14 49 in general is not shape regular In our applications cf chapters 12 and 13 T contains a significant number of strongly deteriorated triangles that have very smal
28. corresponding VertexCL EdgeCL FaceCL and TetraCL objects Note that for a single simplex maybe multiple such num bers have to be stored namely one for each index or in other words one for each finite element type For an extended finite element space a call to UpdateXNumbering augments the usual numbering also called base numbering by a numbering for the new degrees of freedom induced by extension of the finite element space These numbers are not stored in the UnknownHandleCL objects but in a seperate ExtendedIdxCL object It contains a vector xidx MNI where the entry xidx j either stores the number of the extended DoF belonging to the base DoF j or it TO CHAPTER 8 FUNDAMENTAL CONCEPTS AND DATA STRUCTURES contains a flag that the DoF j is not extended Note that UpdateXNumbering has to be called each time the interface has moved to account for the changed extended DoFs Vector and matrix descriptions A VecDescCL object contains a vector Data of type VectorCL and a pointer RowIdx to the asso ciated index of type IdxDescCL Calling the member function SetIdx idx sets the pointer and resizes the vector to the right dimension Similarly a MatDescCL object contains a sparse matrix Data and pointers RowIdx and ColIdx to the associated row and column indices respectively A call of the member funtion SetIdx ridx cidx sets the pointers and deletes all matrix entries The right dimension of the matrix are set later by Sparse
29. defect correction l decoupling of u p c l decoupling of u p o s l decoupling of u p oc 8 Table 1 1 Overview of numerical methods in DROPS 14 Part I Numerical methods 17 In part I of this guide we will explain the numerical methods implemented in DROPS for the numerical simulation of the one and two phase flow problems described in section 1 1 The methods for spatial and time discretization and for the iterative solution of the resulting nonlinear and linear large sparse systems as summarized in table 1 1 are treated In chapter 2 we discuss the multilevel family of tetrahedral nested meshes that is used in all our simulations In the chapters 3 7 we explain the methods used for the models 1 5 given in section 1 1 The presentation follows the row wise ordering of the methods in table 1 1 18 Chapter 2 Hierarchy of tetrahedral grids We outline the basic ideas of the multilevel grid hierarchy on which our finite element dis cretization method is based We only consider multilevel tetrahedral meshes based on red green refinement strategies cf 20 22 27 The idea of a multilevel refinement and coarsening strategy was introduced in 22 and further developed in 24 27 29 63 64 97 This grid re finement technique is used in UG 95 for an overview we refer to 23 25 Similar techniques are used in several other packages We first introduce a few basic notions We assum
30. definite If it is combined with a Krylov subspace method for nonsymmetric problems for example GMRES this SPD property is not required 16 CHAPTER 8 FUNDAMENTAL CONCEPTS AND DATA STRUCTURES 8 8 Input and output In this section we describe input and output interfaces for different types of data Numerical data Vectors and sparse matrices can be saved to and restored from files by using the input and output stream operators gt gt and lt lt implemented for VectorCL and MatrixCL objects The matrix format is MATLAB compatible which is very useful for computing condition numbers or the spectrum of a matrix Geometrical data The initial triangulation 79 can be read from a mesh file generated with the mesh generator GAM BIT 44 To construct the corresponding multilevel triangulation a ReadMeshBuilderCL object containing the mesh file name is passed to the constructor of the MultiGridCL object Here the concept of the MGBuilderCL class is applied cf Section 8 1 from which ReadMeshBuilderCL is derived Other input file formats can be implemented by adding further ancestors of MGBuilderCL For the input and output of a hierarchy of triangulations M 79 77 we use a self defined file format For saving a MultiGridCL object representing a multilevel triangulation we use a software technique called serialization For this reason the class representing this task is called MGSerializationCL The deserialization is done by the cla
31. hyper graphs can often model the communication in a better way cf 40 Thus the load balancing algorithm can be improved by a better node weight function a v cf Definition 12 and by the use of hyper graphs to compute a graph partitioning Reparametrization of the level set function The currently implemented strategy to reparametrize the level set function by performing the update phase in the fast marching algorithm on a single processors is not efficient for a large number of processors Amdahl s law A better parallel variant of the FMM has to be developed and implemented XFEM The present parallel version of DROPS cannot handle extended finite element spaces As is known from the sequential version these spaces result in much smaller discretization errors as the standard finite element spaces Therefore it is important to be able to use the XFEM technique in the parallel code too 88 CHAPTER 9 PARALLELIZATION Part III Examples of implementations in DROPS 91 In this part we present examples of implementations of a few one and two phase flow problems in DROPS 92 Chapter 10 How to get started In this chapter we describe how the DROPS package can be obtained and installed 10 1 How to obtain the code For access to the DROPS CVS server it may be necessary to register at TIM Tivoli Identity Manager at the Computing Center RWTH Aachen University see for details http www rz rwth aachen de Then the
32. mg ZeroFlowCL C DROPS StokesBndDataCL num bnd bc bnd fun DROPS P1X_FE C XFEMStab The function Strategy is now called to solve the problem Strategy prob adap 12 2 4 The function Strategy The Strategy function has the prototype template lt class Coeff gt void Strategy InstatNavierStokes2PhaseP2P1CL lt Coeff gt amp Stokes AdapTriangCL amp adap The triangulation from Stokes is obtained using the routine GetMG 106 CHAPTER 12 NAVIER STOKES EQUATIONS FOR A TWO PHASE FLOW MultiGridCL amp MG Stokes GetMG An object 1set of LevelsetP2CL is created to describe the levelset equation sigma Stokes GetCoeff SurfTens LevelsetP2CL lset MG amp sigmaf grad sigma 0 C lset theta C lset_SD 1 C lset iter C 1set tol C CurvDiff The third parameter in the above constructor of lset is set as the null pointer because the surface tension is constant Otherwise one should use a pointer to the function which describes the gradient of the surface tension Then we create the numberings for the levelset velocity and pressure using the member functions CreateNumbering of lset and Stokes with the corresponding index description pointers lidx vidx and pidx IdxDescCL lidx amp lset idx IdxDescCL vidx amp Stokes vel idx IdxDescCL pidx amp Stokes pr idx lset CreateNumbering MG GetLastLevel lidx Stokes CreateNumberingVel MG GetLastLevel vidx Stokes CreateNumberingPr MG GetLastL
33. n n For each face of the cube the type of boundary condition is set with DirBC which indicates the inhomogeneous Dirichlet boundary condition and the boundary value for the velocity is contained in the array bnd_fun const DROPS BndCondT bc 6 DROPS DirBC DROPS DirBC DROPS DirBC DROPS DirBC DROPS DirBC DROPS DirBC const DROPS StokesVelBndDataCL bnd_val_fun bnd_fun 6 i amp uD amp uD amp uD amp uD amp uD amp uD We define an object prob of the class NavierStokesP2P1CL for the Navier Stokes problem typedef DROPS NavierStokesP2P1CL lt StokesCoeffCL gt NavierStokesCL NavierStokesCL prob brick StokesCoeffCL DROPS StokesBndDataCL 6 bc bnd_fun The problem is solved by calling the function Strategy Strategy prob fp_tol fp_maxiter deco_red oseen_maxiter theta num_timestep time_begin time_end shell_width c_level f_level One can now access the computed solutions via the member functions GetVelSolution and GetPrSolution of prob 11 2 4 The function Strategy The Strategy function has the following prototype template lt class Coeff gt void Strategy DROPS NavierStokesP2P1CL lt Coeff gt amp NS double fp tol int fp maxiter double deco red int oseen maxiter double theta DROPS Uint num timestep double time begin double time end double shell width int c level int f level First we refine the grid in the above mentioned subdomains further and create the triangulat
34. numerical data are coupled The geometrical and numerical data structures are described in Sections 8 1 and 8 2 respectively 8 1 Geometrical objects multilevel triangulation and simplices In this section we discuss the data structures that represent geometrical objects such as ver tices edges faces tetrahedra the boundary and the multilevel grid The corresponding data structures are called VertexCL EdgeCL FaceCL TetraCL BoundaryCL and MultiGridCL re spectively Note that all C classes in DROPS have a suffix CL to distinguish data type identifiers from object identifiers Boundary and boundary segments We assume that the boundary X OQ is partitioned into elementary boundary segments j 0 Ny 1 Note that we use a C style numbering starting with zero To give an example if Q is a cube then X can be partitioned into Ny 6 boundary segments X5 U5 cf Figure 8 2 Each boundary segment is represented by a BndSegCL object Up to now DROPS can only handle boundaries which are piecewise planar The class BoundaryCL contains an array of all BndSegCL objects 8 1 GEOMETRICAL OBJECTS MULTILEVEL TRIANGULATION AND SIMPLICES 67 Simplices In the following we describe the representation of the simplices VertexCL Each vertex V stores its coordinates xy R as a Point3DCL object If V is located on the boundary X it stores a list of BndPointCL objects each containing the index j of the boundary segment X wit
35. of the Oseen solver Concerning the iterative solvers for linear equations Oseen problem and level set equation we distinguish between Schur complement methods only for saddle point problems Krylov subspace methods and multigrid solvers Schur complement methods cf sections 3 5 2 and 3 5 3 e InexactUzawaCL Algorithm 3 28 an a variant of this introduced by Bramble and Pasciak in that is implemented in UzawaCL e SchurSolverCL a variant of algorithm 3 21 3 23 Some of the classes provide template parameters ULPcT SchurPcT to determine the type of the preconditioners QA Qs for the upper left block in the saddle poiint matrix and its Schur complement respectively Krylov subspace methods For the application of a general Krylov subspace method to the saddle point matrix K one can use the class BlockMatrixSolverCL lt SolverT gt where the template parameter SolverT specifies the type of the Krylov solver The following methods are available in DROPS e PMResSPCL preconditioned MINRES solver for the Stokes problem e CGSolverCL PCGSolverCL CG method and preconditioned variant e MResSolverCL PMResSolverCL MINRES method and preconditioned variant e GMResSolverCL GMResRSolverCL GMRES and GMRES Recursive method with left or right preconditioning e BiCGStabSolverCL preconditioned BiCGSTAB method e GCRSolverCL preconditioned GCR method The classes representing preconditioned Krylov subspace methods have a templ
36. of the corresponding quadrature rule Arithmetic operations such as for GridFunctionCL objects are defined point wise In the same way functions can be applied to GridFunctionCL objects using the member function apply Due to inheritence all this functionality is also provided for the derived LocalP CL and Quad CL classes This is very useful when treating complex integrands like u Voj vj Several variants of assign member functions enable the initialization of the LocalP CL and Quad CL objects Additionally LocalP CL objects can be evaluated in an arbitrary point x T given by its barycentric coordinates The Quad CL objects have a member function quad which applies the quadrature rule and returns the result of the numerical integration 72 CHAPTER 8 FUNDAMENTAL CONCEPTS AND DATA STRUCTURES Local numberings A LocalNumbCL object is initialized with an index description of index J the corresponding boundary data object and a tetrahedron T It collects the numbering of the local degrees of freedom of T according to the index J cf Section 8 3 If a degree of freedom is on a boundary it also provides the associated boundary condition and the number j of the corresponding boundary segment X Up to now LocalNumbCL can only be used for P finite elements Integration over interface patches or parts of a tetrahedron An InterfacePatchCL object is initialized by a tetrahedron T and the level set function
37. of this result is given in 1 In 1 a parallel refinement and coarsening algorithm is presented which is based on an admissible hierarchical decomposition and is suitable for distributed memory machines It uses the DDD package 64 95 for the management of the distributed tetrahedra faces edges and vertices For a given input multilevel triangulation the parallel method produces the same output multilevel triangulation as the serial method presented in 27 28 In this sense the computational part of the algorithm is not changed It is proven that the application of the parallel refinement algorithm to an admissible hierarchical decomposition is well defined and results in an admissible hierarchical decomposition Remark 2 Let T M p be a parent master element From the second result in Lemma 1 and A4 it follows that either all children are masters on the same processor p as T or they are masters on some other processor q In the latter case the element T has a corresponding ghost element on processor q Due to this property in the parallel refinement algorithm we use the strategy e If a parent tetrahedron T has a ghost copy then operations that involve children of T are performed on the processor on which the ghost and the children are stored From condition A4 it follows that a child master element has its parent as ghost or as master on the same processor Therefore we use the strategy e Operations that involve t
38. output h include out ensightOut h include geom builder h include stokes stokes h include num nssolver h include navstokes navstokes h include navstokes integrTime h include lt fstream gt include lt sstream gt include lt string gt The program consists of e A class StokesCoeffCL which contains the coefficients of the PDE and the right hand side function g to be used for the template parameter of the problem class NavierStokesP2 P1CL e The functions uD for the Dirichlet boundary conditions and Null for the initial condition of the velocity 11 22 IMPLEMENTATION 99 e The functions MakeInitialTriangulation and ModifyGridStep for generating the initial triangulation e The main function which describes the Navier Stokes problem and calls the function Strategy to solve it e The function Strategy the most important part in which the Navier Stokes problem defined in the main function is solved In the next sections we present the essential parts of functions main and Strategy in detail 11 2 3 The function main The main function first reads the input parameters into the corresponding variables The domain is partitioned uniformly into n x n x n cubes then each of them is dived into 6 tetrahedra using an object of the BrickBuilderCL class DROPS BrickBuilderCL brick DROPS std_basis lt 3 gt 0 DROPS std_basis lt 3 gt 1 DROPS std_basis lt 3 gt 2 DROPS std_basis lt 3 gt 3 n
39. pp 463 502 2001 69 Paige C C Saunders M A Solution of sparse indefinite systems of linear equations SIAM J Numer Anal 11 1974 197 209 70 ParaView parallel visualization application http www paraview org 71 S B PILLAPAKKAM P SINGH A level set method for computing solutions to viscoelastic two phase flow J Comput Phys 174 pp 552 578 2001 12 ParMETIS parallel graph partitioning package http glaros dtc umn edu gkhome metis parmeti 73 Powell C Silvester D Black box preconditioning for mixed formulation of self adjoint elliptic PDEs In Challenges in Scientific Computing CISC 2002 ed E B nsch Lecture Notes in Computational Science and Engineering Vol 35 268 285 2003 74 R RANNACHER Finite Element Methods for the Incompressible Navier Stokes Equations in Fundamental directions in mathematical fluid mechanics G P Galdi et al eds pp 191 293 Birkh user Basel 2000 75 R RANNACHER On the numerical solution of the incompressible Navier Stokes equations ZAMM 73 pp 203 216 1993 V C Boffy and H Neunzert eds pp 34 53 Teubner Stuttgart 1998 76 H G Roos M STYNES L TOBISKA Numerical methods for singularly perturbed differ ential equations convection diffusion and flow problems Springer ser in comp math Vol 24 Springer Berlin Heidelberg 1996 77 M RUDMAN Volume Tracking Methods for Interfacial Flow Calculations Int J Num Me
40. pp given by an P5 FE VecDescCL object It extracts the LocalP2CL object corresponding to yp decides whether T N T Z and provides information about the sign 0 of each degree of freedom The member function ComputeForChild i computes the planar interface patch rr TRN T for the ith regular child T K T i 0 7 Tr is represented by the coordinates of its vertices which are given in terms of barycentric coordinates with respect to the parent T cf Fig 4 1 and Fig 4 2 Note that for the computation of the patches the regular refinement of T is not really constructed in the sense that geometrical data structures are changed After calling ComputeCutForChild i the member function quad can be used to compute the integral over the subdomain T N Q or T n Q2 where the integrand is an ar bitrary quadratic function f given by a LocalP2CL object The additional member function quadBothParts provides the integrals over the union of the subdomains T N Q1 U TN Q2 8 6 Time discretization and coupling For the one phase Stokes and Navier Stokes problem the one step 0 scheme cf section 3 3 1 is represented by the classes InstatStokesThetaSchemeCL and InstatNavStokesThetaSchemeCL respectively Both classes have a template parameter SolverT for the type of the solver used in each time step The computation of one time step is performed by the member function DoStep For the two phase Stokes and Navier St
41. problem in Qa C H1 Q Trx y VJox V Joy for all x y RF 3 5 ITERATIVE SOLVERS FOR LARGE SPARSE LINEAR SYSTEMS 29 with Jo R Qn the finite element isomorphism between finite element functions in Qa and their nodal values cf 3 4 Note that ker T span e with e 1 1 1 T We define Q et RX by 2 M 6T if B lt h Q Se 3 32 BK M BT if h g Note that for 8 0 which corresponds to the discrete stationary Stokes equation we get Qs M In 34 it is shown that Qs is uniformly in and spectrally equivalent to the Schur complement S Preconditioned MINRES method We consider a preconditioned minimal residual PMINRES method for solving the discretized Stokes problem 3 20 This class of methods has been analyzed in 73 78 85 We consider a symmetric positive definite preconditioner k Qa 0 K 3 33 eE 3 33 for K Define the norm v z Kv v 8 for ve REN Given a starting vector v with cor responding error e v v9 then in the preconditioned MINRES method one computes the vector v v span K Ke HEP K K e which minimizes the preconditioned residual K K v v over this subspace An efficient implementation of this method can be derived using the Lanczos algorithm and Givens rotations For such an implementation we refer to the literature e g 69 51 In an efficient implementation of this method one needs per iteration one evalu
42. sharp w r t the order of convergence for hr 0 Using the second Strang lemma cf 4 this approximation error induces an error of the same order of magnitude in the H norm in the velocity approximation and thus relatively large spurious velocities can occur In 3 a modified surface tension force discretization with better approximation quality is presented We briefly explain this method For this we have to introduce some further notation Let n be the unit normal on T outward pointing from Q1 Since I is planar piecewise triangular this normal is piecewise constant and not defined at the edges of the surface triangulation We define the orthogonal projection Py P z I n z n 2 for zx T x not on an edge Recall that the discrete level set function j is piecewise quadratic Define z Von 2 I Von 2 For the discrete surface tension force as in 4 14 we have due to Vr g P Vg for smooth functions g the representation P z I n z n 2 x T x not on an edge 3 fr Vn DL P x e Vr vn i ds 4 16 36 CHAPTER 4 NAVIER STOKES EQUATIONS FOR TWO PHASE FLOW NS2 with e the i th basis vector in IR and vn i the i th component of vp The modified discrete surface tension force is given by 3 fr vn p P z e Vr va i ds 4 17 The implementation of this functional requires only a minor modification if the implementation of the one in 4 16 is available
43. t lt lt t lt lt maxGradPhi lt lt t lt lt Volume lt lt t lt lt bary drop lt lt t lt lt min drop lt lt t lt lt max drop lt lt std endl std cerr lt lt rel Volume lt lt lset GetVolume Vol lt lt std endl if C VolCorr double dphi lset AdjustVolume Vol 1e 9 Std cerr lt lt volume correction is lt lt dphi lt lt std endl lset Phi Data dphi std cerr lt lt new rel Volume lt lt lset GetVolume Vol lt lt std endl if C RepFreq amp amp stepAC RepFreq 0 reparam levelset function lset ReparamFastMarching C RepMethod adap UpdateTriang Stokes lset amp c if adap WasModified 1 cpl Update c Update std cerr lt lt rel Volume lt lt lset GetVolume Vol lt lt std endl if C VolCorr 1 double dphi lset AdjustVolume Vol 1e 9 Std cerr lt lt volume correction is lt lt dphi lt lt std endl lset Phi Data dphi std cerr lt lt new rel Volume lt lt lset GetVolume Vol lt lt std endl ensight putGeom datgeo step C dt ensight putScalar datpr Stokes GetPrSolution step C dt ensight putVector datvec Stokes GetVelSolution step C dt ensight putScalar datscl lset GetSolution step C dt ensight putScalar datc c GetSolution step C dt ensight putScalar datct c GetSolution c ct step C dt ensight Commit 13 3 Result
44. tool the tetrahedra and numer ical data are rearranged among the processors This phase is called data migration To obtain again an admissible hierarchical decomposition after the migration phase we have to ensure that the properties A1 A5 hold In particular all children of a common parent have to stay to gether as masters on a single processor cf Lemma 1 Thus in the following we give a definition for a reduced dual graph where the children of a common parent are represented by a single multi node For this purpose we introduce a map J J i P U Gk gt U Gk k 0 k 0 from a tetrahedron T G to its parent tetrahedron P T G _1 k 1 J with the convention P T T for all T Go For T 7 we define the corresponding equivalence class Tp S eT PS P r 84 CHAPTER 9 PARALLELIZATION Figure 9 4 Reduced dual graph for 2D triangulation Definition 12 Reduced dual graph For a triangulation 7 let G T V E o 3 be the corresponding weighted dual graph The reduced dual graph G 7 V E a B is given by the reduced node set V T p TET inducing the reduced edge set E vi vh Jur Evl v E vy vv e EJN Q v v eV The weight functions a 3 are given by e En vj x v5 Figure 9 4 shows the reduced dual graph corresponding to the dual graph given in Figure 9 3 The tetrahedra forming a multi node are surrounded by a bold frame Note that the dual graph G T
45. wa 0 Ot divu 0 in the unit cube Q 0 1 for t 0 1 with u 0 001 This rather small viscosity coefficient is used to obtain a more interesting compared to u 1 flow field cf Figure 11 2 The Dirichlet boundary conditions are taken as u 1 0 0 7 on the face z 0 and u 0 0 0 7 on the other faces The initial condition is u 0 0 0 7 for t 0 The following components are used e For the triangulation we partition Q uniformly into 16 x 16 x 16 cubes and then each of them is subdivided into six tetrahedra Then the subdomains 0 0 2 x 0 0 2 x 0 0 2 and 0 8 1 x 0 0 2 x 0 0 2 are refined once more e P P finite element pair for velocity and pressure The discretization is implemented in the routines SetupInstatSystem SetupInstatRhs of the class StokesP2P1CL file stokes stokes h and SetupNonlinear of the class NavierStokesP2P1CL file navstokes nav stokes h e For time discretization and linearization the 0 scheme with 0 1 and the fixed point de fect correction are used cf sections 3 3 1 and 3 4 These are implemented in the classes In statNavStokesThetaSchemeCL file navstokes integrTime h and AdaptFixedPtDefect CorrCL file num nssolver h respectively e The Oseen problems are solved by a preconditioned block GCR solver cf section 3 5 3 In DROPS it is realized by a combination of the classes GCRSolverCL file num solver h and BlockMatrixSolverCL file num stokessolver h The acc
46. 1 is a reformulation of the 0 scheme applied to the system of ODEs in 4 34 The latter is A stable and first order accurate for 0 0 1 and second order accurate for 0 3 Remark 7 For the case 0 1 the scheme takes a much simpler form In particular the sequences for Z and w are not needed The resulting method is as follows R ES 70 A Given u determine for n gt 0 2ncl gt n 1 y 4 gt n 1 gt 2n F1 E yu N at yan g m fr 4 B p gt n 1 ur M PAS Bu t 0 gt n 1 gt n n 1 p E u AS E d HH x 4 3 2 The generalized 0 scheme for a time dependent B In this section we allow both M and B to be time dependent We use the same notation as in the previous section For the derivation of a generalized 0 scheme we apply an approach that is 4 3 TIME INTEGRATION 4T standard in the field of differential algebraic equations cf We introduce gt M 9 G g fr BT p Fi y y C y E u F y H u F2 y p 0 B F y The Navier Stokes equations coupled with the level set equation can be written as dy CT F y with initial condition for to and to For the numerical treatment of this DAE system we introduce the variable q ay Y and thus we obtain the coupled DAE system dy E le a 0 d 4 43 C y q F y 0 The 0 scheme with 0 Z 0 applied to this takes the form y n At 6d qi
47. 2DSolOutCL for visualization of geometry and numerical solu tion in 3D or on a 2D cut plane respectively with TecPlot 89 e MapleMGOutCL MapleSolOutCL for visualization of geometry and numerical solution with Maple Chapter 9 Parallelization In this chapter we consider the main concepts underlying the parallelization of DROPSfor dis tributed memory machines by means of a message passing interface MPI Shared memory parallelization by means of OpenMP 66 has also been applied to some parts of DROPS cf 14 for a description of the parallelized routines and some benchmark computations Both paral lelization concepts can be combined when using multicore processors which are connected by a high speed network For the parallelization of DROPS we pursue such a hybrid parallelization approach due to the growing importance of multicore architectures In Section 9 1 we present a data distribution format for the geometrical data and based on this we also derive a distribution format for the numerical data In Definition 7 below the geo metrical data distribution format will be made mathematically precise by a formal specification of a so called admissible hierarchical decomposition lhis data distribution format is such that the following holds 1 Let T Gk be an element from the hierarchical surplus on level k cf Definition 6 Then T is stored on one processor say p as a so called master element In certain cases explained belo
48. 58 Two phase flow model 1 4 combined with btn divr Sur SKu n DrArS 1 6 The derivative 0 S stands for the time derivative of S along a normal path The diffusion coefficient Dp can be assumed to be constant on I For the convection diffusion equation in 1 6 no boundary conditions are needed if the interface T is a manifold without boundary 5 NS2 combined with transport of both a dissolved species and a surfactant at the interface NS24 T4 S Forthcoming This model combines the models in 1 5 and 1 6 If c in 1 5 models the concentration of the surfactant in the two phases then in 1 6 one usually includes an additional source term that accounts for the change of the surfactant concentration at the interface S due to ad and desorption cf 31 98 12 CHAPTER 1 INTRODUCTION 1 2 Initial and boundary conditions In this section we describe the initial and boundary conditions that can be used in the models 1 5 to make the problem well posed For the NS1 model one needs suitable initial and boundary conditions only for the velocity u The initial condition is u x t uo r with a given function ug which usually comes from the underlying physical problem For the boundary conditions we distinguish between essential and natural boundary conditions Let ON be subdivided into two parts Q O p U ONN with Np N 0 y 0 We use essential boundary conditions on Np that are of Dirichlet type In application
49. DROPS Package User s guide March 19 2008 Contents 1 Introduction 1 1 One and two phase flow models in strong formulation 1 2 Initial and boundary conditions CE Emm nn 1 3 Overview of numerical methods o oaoa I Numerical methods 2 Hierarchy of tetrahedral grids 3 Navier Stokes equations for one phase flow NS1 3 1 Weak formulation 4 4253 24044 8 boa Sr ee do d UA ORO I ESAE ROC 3 2 Spatial discretiZatiOncz us ule Boa Ree neck qe el pedes YO 3 2 1 3 2 2 Galerkin finite element discretization enn Quadrature Ya Wa aa SCHON re S RUBIO EL MS cs 3 5 Lime integration 4 224 Lg ko ee wel EA ae eo ea 3 9 1 3 3 2 3 3 3 Th sheme ahnen te an Ex oe BO ee E EIE RES fractional step 6 scheme 2h Operator splitting les 3 4 Linearization methods lll ln 3 5 Iterative solvers for large sparse linear systems 00000 ae 3 5 1 3 5 2 3 5 3 Iterative solvers for discretized scalar elliptic problems Iterative solvers for discretized Stokes equations Iterative solvers for discretized Oseen equations 4 Navier Stokes equations for two phase flow NS2 4 1 Weak formulations are ko Eck eb sepa ee ea al dedo 4 2 Spatial discretization 2 2 uoa gen Bt eB eh en aSa 4 2 1 4 2 2 4 2 3 4 2 4 4 2 5 4 2 6 4 2 7 Galerkin finite element discretization ee eee Disc
50. In 3 it is shown that under reasonable assumptions on I and op we have the error bound sup 1fe va fex Cv l fr va lt chp 4 18 vn amp eVn Iva lli This bound has a O hr behaviour instead of O hp Numerical experiments in 3 show a rate of convergence that is even somewhat higher than first order in hp In DROPS both discretizations 4 14 and 4 17 of the localized force term fr are available The default choice is the modified functional in 4 17 Remark 3 The case with a variable surface tension coefficient T T x can also be treated in DROPS Using partial integration the surface tension functional in 4 2 takes the form fr v nz Vpr vds i Kn rv Vr v ds 4 19 3 J T Vridr Vrv ds gt Vr idr Vrr vi ds VrT v ds r SUP T For discretization of this term we can easily generalize the approach discussed above for the case of a constant surface tension coefficient The interface I is replaced by its approxmation I and the tangential derivative are approximated using the projection P in the same way as in 4 17 This results in the discrete surface tension functional 3 f v P i V v id fe va er a z ei Vr va ds 3 i 2 P x ei Vr 7 vi i ds ji Py z Vr vn ds Th We assume that 7T T x which is defined for x L has an extension to z Tp 4 2 4 Modified finite element space for the pressure XFEM In this section we discuss t
51. J Comp Phys 114 pp 146 159 1994 M SUSSMAN A S ALMGREN J B BELL PH COLELLA L H HOWELL M L WELCOME An adaptive level set approach for incompressible two phase flows J Comp Phys 148 pp 81 124 1999 Tecplot Tecplot Inc 2D 3D visualization tool http www tecplot com A A JOHNSON T E TEZDUYAR Mesh update strategies in parallel finite element compu tations of flow problems with moving boundaries and interfaces Comput Methods Appl Mech Engrg 119 pp 73 94 1994 A K TORNBERG Interface tracking methods with application to multiphase flows Doc toral Thesis Royal Institute of Technology Department of Numerical Analysis and Com puting Science Stockholm 2000 A K TORNBERG B ENGQUIST A finite element based level set method for multiphase flow applications Comput Visual Sci 3 pp 93 101 2000 G TRYGGVASON B BUNNER A ESMAEELI D JURIC N AL RAWAHI W TAUBER J HAN S NAs Y J JAN A front tracking method for the computations of multiphase flow J Comp Phys 169 pp 708 759 2001 S TUREK Efficient solvers for incompressible flow problems An algorithmic approach in view of computational aspects LNCSE Vol 6 Springer Berlin Heidelberg 1999 UG http cox iwr uni heidelberg de ug S O UNVERDI G TRYGGVASON A front tracking method for viscous incompressible multi fluid flows J Comp Phys 100 pp 25 37 1992 S ZHANG Successive subdivisions of tetrahedra
52. MM cf section 4 2 6 is a delicate task because information from neighboring tetrahedra is needed to update degrees of freedom This information is used in each update step of the FMM and may be stored on a neighbor processor This implies a large number of small communications and thus can cause bad parallel efficiency Due to the inherent sequential structure of this algorithm there is no effective way of parallelizing this algorithm with only minor changes We have implemented in the current version a with respect to parallel efficiency suboptimal strategy in which the update phase of the FMM is realized on a single processor It turns out that this is still acceptable for the problem sizes and numbers of processors that we have considered so far For the decoupling in each time step between the level set and the Navier Stokes subprob lems in a two phase flow problem the sequential methods could be parallelized with only minor changes Outlook on parallelization The following concerning the parallelization of DROPS are planned e Hybrid parallelization In order to use clusters of multicore architectures a hybrid parallel version of DROPS will be developed This will be done by combining OpenMP and MPI to get better parallel efficiency on these high performance computers The planned shared memory parallelization will be based on the work that has been done by the Center for Computation and Communication at RWTH 14 e Improving paral
53. MatBuilderCL cf Section 8 2 8 4 Problem classes There are several problem classes in DROPS representing different types of partial differential equations For example we have problem classes for the Poisson Stokes and Navier Stokes problem one phase the level set equation and the two phase Stokes and Navier Stokes problem For example the class for the two phase Stokes problem is called InstatStokes2PhaseP2P1CL All problem classes are derived from a common base class ProblemCL which contains three objects constituting a problem e the domain 2 given by a multilevel triangulation MultiGridCL e the boundary conditions and boundary values given by a BndDataT object e the coefficients and right hand side of the partial differential equation given by a CoeffT object BndDataT and CoeffT are template parameters of the template class ProblemCL as their specific structure may vary among different problem types Their meaning is discussed in the subsequent sections A specific problem class usually contains the index descriptions of the applied finite ele ment types and several matrix and vector descriptions Among the member functions there are CreateNumbering procedures for the indices cf Section 8 3 and different Setup routines to compute the matrices and the right hand side vectors constituting the finite element discretization Boundary data The boundary data are represented by a BndDataCL lt BndValT gt object
54. N P Ps Pi T VA en y A K T A AA NAN T A B Figure 9 1 Ghost elements are required to represent links between parents and their children since pointers across memory boundaries are not allowed for distributed memory machines e Given an admissible hierarchical decomposition one can apply a suitable load balancing and data migration algorithm such that after data migration one still has an admissible hierarchical decomposition We comment on this in Section 9 2 9 1 Data distribution Distribution of geometrical data admissable hierarchical decomposition Let the sequence M To Ty of triangulations be a multilevel triangulation and H Co GJ the corresponding hierarchical decomposition In this section we introduce a par ticular format for the distribution of the tetrahedra in H among processors on a parallel machine We assume that the processors are numbered by 1 P For the set of elements in the hierarchical surplus on level k that are stored on processor p we introduce the notation Gi p T Gy T is stored on processor p and we define Hlp o p 95 p Note that in general H p is not a hierarchical decomposition in the sense of Definition 7 The sequence H H 1 H P 9 1 is called a distributed hierarchical decomposition corresponding to H In general the intersection Gk p N Gr q p q may be nonempty Note that such an overlapping distribution of th
55. OKES EQUATIONS FOR A ONE PHASE FLOW Figure 11 1 Driven Cavity Velocity field at steady state Figure 11 2 Driven Cavity Streamlines at steady state Chapter 12 Navier Stokes equations for a two phase flow 12 1 Introduction We consider a model of the form 1 4 that describes a rising bubble problem namely an n butanol droplet in water that rises due to gravity The domain is given as Q 0 0 008 x 0 0 04 x 0 0 008 At the initial time the droplet 01 is a sphere of radius R 0 001 centered at 0 004 0 002 0 004 The material coefficients are given as pii 0 003281 u2 0 001388 p 845 442 12 1 p2 986 506 T 0 00163 Homogeneous Dirichlet boundary conditions for velocity are used The initial velocity is 0 0 0 7 and the external gravity force is 0 9 81 0 7 The following methods are used e P P X finite element spaces for velocity and pressure The discretization of the two phase flow Navier Stokes equations is done in the class InstatNavierStokes2PhaseP2P1CL with the routines SetupSystemi SetupSystem2 and SetupNonlinear in the files stokes in statstokes2phase h and stokes instatnavstokes2phase h e The spatial discretization of the levelset equation is implemented in the classes Levelset P2CL and InterfacePatchCL file levelset le
56. PcPcT APcPcT Apcpc typedef GMResSolverCL lt APcPcT gt ASolverT ASolverT Asolver Apcpc 500 restart 100 1e 2 relative true typedef SolverAsPreCL lt ASolverT gt APcT APcT Apc Asolver typedef MinCommPreCL BBTPcT BBTPcT BBTpc 0 Stokes B Data Stokes M Data Stokes prM Data typedef BlockPreCL lt APcT BBTPcT false gt OseenPcT OseenPcT oseenpc Apc BBTpc typedef GCRSolverCL lt OseenPcT gt OseenBaseSolverT OseenBaseSolverT oseensolver0 oseenpc truncate C outer_iter C outer iter C outer tol relative false typedef BlockMatrixSolverCL lt OseenBaseSolverT gt OseenSolverT OseenSolverT oseensolver oseensolverO0 typedef AdaptFixedPtDefectCorrCL lt StokesProblemT OseenSolverT gt NSSolverT NSSolverT nssolver Stokes oseensolver C ns_iter C ns_tol C ns_red An object of the class LinThetaScheme2PhaseCL is created for time integration typedef LinThetaScheme2PhaseCL lt StokesProblemT NSSolverT gt CouplingT CouplingT cpl Stokes lset nssolver C nonlinear true C cpl_stab The routine SetMatrixA is used to set the actual matrix A to the BBTPcT preconditioner BBTpc SetMatrixA amp cpl GetUpperLeftBlock We use the routine SetTimeStep to set the time and the time integration procedure is as follows cpl SetTimeStep C dt for int step 1 step lt C num_steps step 1 cpl DoStep C cpl iter time C dt if C VolCorr 1 double dphi lset AdjustVolume
57. TAB GCR 3 5 ITERATIVE SOLVERS FOR LARGE SPARSE LINEAR SYSTEMS 27 For a description of these methods we refer to the literature e g 54 79 These methods which are easy to implement are available in DROPS However often their rate of convergence is very low A much more efficient iterative solver can be obtained by using basic iterative methods like a damped Jacobi or Gauss Seidel method in a multilevel approach Since in our discretization technique we have a hierarchy of nested grids and corresponding finite element spaces available such a multigrid approach can be used A multigrid iterative solver with a damped Jacobi or Gauss Seidel smoother and canonical intergrid transfer operators induced by the nestedness of the finite element spaces has been implemented in DROPS This method is a very efficient solver for scalar diffusion and convection diffusion problems 3 5 2 Iterative solvers for discretized Stokes equations The discrete Stokes problem has a matrix vector representation of the form K E JG R Anatom 3 20 where the parameter is proportional to 1 At The matrix is symmetric positive definite The matrix K is symmetric and strongly indefinite and has a saddle point structure For this type of linear systems we consider the following methods e inexact Uzawa methods e preconditioned MINRES e multigrid method We briefly address these three classes of methods Inexact Uzawa method The Schur com
58. a 1 g 42 SPATIAL DISCRETIZATION A1 Another approach called fast marching method FMM is based on a greedy grid traversal technique combined with a function for computing approximate distances Such an FMM has been implemented in DROPS We outline this method It consists of two phases an initialization phase where the vertices that are directly adjacent to I are assigned new values and a second far field phase in which all the vertices further away are updated Let I be the discrete interface that has been constructed as described in section 4 2 2 and Ty the uniform refinement of Tp that is used in this construction Let Ty T T meas3 T T4 gt 0 be the set of tetrahedra that have a nonzero intersection with the interface The discrete interface Ta consists of planar segments and can be represented as Tl Ure Tr with Ir D4nT For a vertex v the set of tetrahedra T 7 that have v as a vertex is denoted by T v For T Tj the set of vertices of T is denoted by V T Initialization phase In the initialization phase an approximate disctance function for the vertices v Vp v Y T T Tp is determined For a given T Tr let rr I NT be the planar segment that is part of I cf Fig 4 2 and let Qi Qm m 3 or 4 be the vertices of this segment Let W CR be the plane that contains Tr and Py the orthogonal projection on W For T Tp and v V T we define iu l
59. a 99 WS Results sinas ana DR ero us ton tee Tees ek pee UR ale ahd o p ot tee Gore ek 101 12 Navier Stokes equations for a two phase flow 103 1221 IHtr duction s 2 a a cabs a mox oO eee Nr en SA een en eg 103 12 2 Implementation 4 2 22 25 A aa sr He ee Dae Be ee s 104 12 3 1 Input parameters s 28 ego om RR ek ee ee eke A 104 12 2 2 Structure of the program mn nennen 104 12 2 3 The function main 2 222 rs 105 12 2 4 Th function Strategy vii os 54 reo aa Xe Rb RI OR nr Warn 105 12 3 Resulbss Ass 8 3 thee expo ete Lo ee ru de send oid unb e sop OS ee uis 108 13 A two phase flow problem with mass transport 111 13 12 Introductions Ara hut uds die nist Babee ped ho Seduta d 111 13 2 Implementation ooa ee 112 13 37 esu Sea dos VoM ck Om em zs e entes SU Sd Mr UAR E UM Ep EE a ca s Ust ua 113 14 Appendix 115 LAST Parameter fileant 8 3 28 3 ete CROSS afud eot CU A 115 CONTENTS Chapter 1 Introduction The development of the software package DROPS started in the interdisciplinary Collaborative Research Center SFB 540 Model based Experimental Analysis of Kinetic Phenomena in Fluid Multi phase Reactive Systems 84 The goal of the DROPS related research activities is two fold on the one hand we want to develop and improve numerical methods for the simulation of two phase flow problems and on the other hand the aim is to provide accurate and reliable numerical simulations of realistic
60. ab E poe Tecplot ASCII e5 Geomview 7O I b 4 p Geometry Numerics Figure 8 1 Overview of modules and structure of DROPS 65 66 CHAPTER 8 FUNDAMENTAL CONCEPTS AND DATA STRUCTURES x Figure 8 2 A cube and its 6 boundary segments X5 U5 The vertical structure of the figure distinguishes between input and output routines data structures and algorithms The different methods related to input and output are described in Section 8 8 In part I we treated many numerical methods from a numerical analysis point of view In this chapter we discuss implementation issues of these algorithms in particular their relation with certain data structures This corresponds to the object oriented perspective of C classes where data structures as data members and functionality as member functions are combined with each other The horizontal structure in Figure 8 1 shows a classification of the different modules based on the categories geometry and numerics emphasizing the fact that we tried to decouple geo metrical data such as the grid from numerical data such as vectors and matrices Some tasks however require geometrical as well as numerical information and are therefore located in the middle column For example the discretization routines for setting up stiffness matrices where in a loop over all tetrahedra the corresponding matrix entries are determined In these routines the geometrical and
61. al data all units are SI Mat DensDrop 845 442 n Butanol Wasser ViscDrop 3 281e 3 DensFluid 986 506 ViscFluid 1 388e 3 SmoothZone 1e 4 SurfTension 1 63e 3 Y mass transport parameters all units are SI MassTransp ConcDrop 5 ConcFluid 0 1 DiffusityDrop 2 3e 6 DiffusityFluid 5 76e 6 EquilibJump 0 75 cFluid EquilibJump cDrop in the equilibrium Tol 1e 10 relative Iter 500 Y experimental conditions Exp RadDrop 1e 3 1e 3 1e 3 PosDrop 4e 3 2e 3 4e 3 Gravity 0 9 81 0 FlowDir 1 flow in y direction InflowVel 0 RadInlet 4e 3 3 5e 3 for old meas cell Y miscellaneous InitialCond 0 O zero 1 2 flow with without drop 3 read from file InitialFile initial brick_adapO MeshFile 8e 3x40e 3x8e 304x20x4 EnsightCase butanol_diam2mm EnsightDir ensight XFEMStab 1 0 Bibliography mm 10 11 12 13 S Gross A REUSKEN Parallel multilevel tetrahedral grid refinement SIAM J Sci Comp 26 No 4 pp 1261 1288 2005 S Gross V REICHELT A REUSKEN A Finite Element based Level Set Method for Two phase Incompressible Flows Comp Vis Sci 9 No 4 pp 239 257 2006 S Gross A REUSKEN Finite element discretization error analysis of a surface tension force in two phase incompressible flows SIAM J Numer Anal 45 No 4 pp 1679 1700 2007 S GROSS A REUSKEN An extended pressure finite element space for two phase incom
62. an curvature flow and related topics In Frontiers in Numerical Analysis Durham 2002 J F Blowey A W Craig T Shardlow eds pp 63 108 Springer 2003 K D DEVINE E G BOMAN R T HEAPHY R H BISSELING U V CATALYUREK Par allel hypergraph partitioning for scientific computing Parallel and Distributed Processing Symposium 2006 IPDPS 2006 20th International 10 pp 2006 G DZIUK An algorithm for evolutionary surfaces Numer Math 58 pp 603 611 1991 H ELMAN V E HowLE J SHADID R SHUTTLEWORTH R TUMINARO Block pre conditioners based on approximate commutators SIAM J Sci Comp 27 pp 1651 1055 2005 Ensight CEI Inc 3D visualization tool http www ensight com GAMBIT ANSYS Inc mesh generation tool http www fluent com S GANESAN G MATTHIES L TOBISKA On spurious velocities in incompressible flow problems with interfaces Preprint Department of Mathematics University of Magdeburg 2005 S GANESAN L TOBISKA Finite element simulation of a droplet impinging a horizontal surface Proceedings of ALGORITMY 2005 pp 1 11 2005 Geomview 3D viewing tool http www geomview org 120 48 49 50 5l 52 53 54 55 56 57 58 59 60 61 62 63 64 BIBLIOGRAPHY I GINZBURG G WITTUM Two phase flows on interface refined grids modeled with VOF staggered finite volumes and spline interpolants J Comput Phys 166 pp 302
63. and multigrid methods on tetrahedral meshes Houston J Math 21 pp 541 556 1995 H ZHOU C CRISTINI J LOWENGRUB C W Macosko Numerical simulation of de formable drops with soluble surfactant Pair interactions and coalescence in shear flow Preprint 1898 Institute of Applied Mathematics University of Minnesota 2002
64. ate parameter PcT designating the type of the preconditioner Multigrid method The MGSolverBaseCL represents a multigrid solver V cycle with a fixed number of smoothing steps There are two template parameters SmootherT and SolverT which control the type of the smoother and the coarse grid solver respectively The multigrid method requires a hierarchy of linear systems Axe b Ellen 8 7 ITERATIVE SOLVERS AND PRECONDITIONERS 75 and prolongations P and restrictions Rg PT to transfer information between a finer and a coarser level For each level the corresponding system and prolongation matrices Ap P and right hand side vector bg are stored in a MGLevelDataCL object The hierarchy of matrices and vectors is represented by the data structure MGDataCL which is simply a list of MGLevelDataCL objects For the smoothers we implemented methods that are suitable for scalar convection diffusion problems such as damped Jacobi and Gauss Seidel and methods that can be used for Stokes and Oseen equations for example Vanka and Braess Sarazin smoothers Preconditioners The DROPS solver toolbox comprises the preconditioner classes given in the following lists For a discussion of some of these preconditioners we refer to sections 3 5 2 and 3 5 3 Iterative method as preconditioners e JACPcCL one step of the Jacobi preconditioner e GSPcCL SGSPcCL one step of the Gauss Seidel or symmetric Gauss Seidel preconditioner e SSORPcCL MultiSSORP
65. ation of QS one eval uation of Qu and one matrix vector product with K From a convergence analysis of this method it follows that we have fast convergence if Q4 and Qs are good preconditioners of A and S respectively For these preconditioners one can take the ones discussed above A detailed discussion of this method and a comparison with the inexact Uzawa method are given in 10 Multigrid method Multigrid methods can be applied directly to the generalized Stokes problem 3 20 These so called coupled multigrid methods are based on a combination of a smoother that is applied to the saddle point system 3 20 and coarse grid corrections that are obtained from discretizations of the form 3 20 on coarser grids For a treatment of these multigrid methods we refer to the literature eg 33 53 In DROPS two types of smoothers are available namely a Braess Sarazin method and a Vanka smoother An explanation of these smoothers and of the corresponding multigrid solver is given in 5 In that paper results of numerical methods are given in which this multigrid method is compared to the preconditioned MINRES and the inexact Uzawa method These results show that for the class of generalized discrete Stokes equations as in 3 20 all three methods are robust with respect to variation in the parameters h and 8 and often these three methods have a comparable efficiency 3 5 3 Iterative solvers for discretized Oseen equations The discrete Oseen prob
66. bble for example the multilevel triangulation M will change as the refinement zone is moving upwards following the bubble geometry Hence the distributed hierarchical decomposition H and the numerical data have to be redistributed from time to time to ensure a reasonable balance of the computational load Otherwise the situation may occur that almost all unknowns are stored on one processor say p while the others only have to solve problems of small size On the one hand this leads to an inefficient usage of the overall memory On the other hand runtime scalability severely decreases since all processors have to wait at synchronization points such as MPI AllReduce until processor p has finished its work The challenge of the so called load balancing is to find a mapping m T gt 1 P describing the distribution of the tetrahedra among the processors such that a all processors have almost the same number of tetrahedra and b the corresponding processor boundary is as small as possible This problem statement is equivalent to a graph partitioning problem which will be stated in Definition 11 For this reason m is also called a partitioning of T We now introduce the notion of a weighted dual graph Definition 10 Weighted dual graph For a triangulation 7 the corresponding dual graph G T V E is given by the node set V 7 and the edge set E C 7 x 7 where T1 T3 E iff the tetrahedra T1 75 share a common face By intro
67. cCL one or multiple steps of the SSOR preconditioner e MGPreCL fixed number of multigrid V cycles with SSOR smoothing e DummyPcCL no preconditioning For the Jacobi and Gauss Seidel preconditioners there exist variants which can be used as smoother for the multigrid solver The wrapper class SolverAsPreCL enables the use of a solver object as a preconditioner That means that the Apply member function of the wrapper class calls the Solve member function of the solver class with initial guess zero This mechanism is used in Exam ple 8 1 in the definition of the preconditioner for the upper left block ULPc which wraps the solver object ULsolver Schur complement preconditioners Qs e ISPreCL the Schur complement preconditioner 3 32 where M and T are replaced by one step of the SSOR preconditioner applied to the corresponding pressure matrices e ISNonlinearPreCL the same Schur complement preconditioner but with M and n replaced by some iterations of a Krylov subspace method which can be chosen by means of a template argument e ISBBTPreCL variant of the Schur complement preconditioner 3 36 The DiagBlockPreCL is used in combination with BlockMatrixSolverCL solvers It com bines a preconditioner Q4 for the upper left block with a preconditioner Qs for the Schur complement yielding the diagonal block preconditioner K as in 3 33 If this preconditioner is used in combination with MINRES it must be symmetric positive
68. ce 59 II Implementation in DROPS 61 8 Fundamental concepts and data structures 65 8 1 Geometrical objects multilevel triangulation and simplices 66 8 2 Numerical objects vectors and sparse matrices 2 oo a 68 8 3 The connnection between grid and unknowns indices 69 8 4 Problem classes ud esie aa e tnb Ir eyed pore Are at Gl 70 8 5 Tools for spatial discretization 2 2 so o a a 71 8 6 Time discretization and coupling 0 000 eee eee 72 8 7 Iterative solvers and preconditioners 2 2m mm mn 73 8 8 Input and output 23 wes Se ae A OE a RA Rd d 76 9 Parallelization 77 9 L Datardistribution m casna 2 2 8 Ros oom ns on in a Pom ig 78 9 2 Distribution of workload es 82 9 3 Current status and outlook on on nn 85 III Examples of implementations in DROPS 89 10 How to get started 93 10 1 How to obtain the code rn 93 10 2 Compilation of DROPS lle s 94 10 3 Building the bimmary files top sono ue ead Rx Wet RR RR Rd em D Res 95 10 4 Running DROPS osos e RR URSUS x en 95 11 Navier Stokes equations for a one phase flow 97 11 1 1ntrod ction s 2 2 2 ma ne Po ae Ses Oe 97 11 2 Implementation 452 Egan Mann sl Drei v 98 11 2 1 Input parameters oe pte 5 6 ee ee son REO ideeller 98 11 2 2 Structure ofthe program aaa 98 CONTENTS 5 1 1 2 3 Phe function main 56 RlIeibU uo 6 8 ite ode dik ek amp ee 2s 99 11 2 4 The function Strategy
69. ce function is discretized with continuous piecewise quadratic finite elements on the tetrahedral triangulation 7 The resulting piecewise quadratic finite element approximation of on T is denoted by n i x t We assume a given fixed time t We introduce one further regular refinement of 7p resulting in 7 Let I amp be the continuous piecewise linear function on 7 which interpolates at all vertices of all tetrahedra in 77 The approximation of the interface T is defined by Ty z Q I s zx 20 4 10 34 CHAPTER 4 NAVIER STOKES EQUATIONS FOR TWO PHASE FLOW NS2 Figure 4 1 Construction of approximate interface for 2D case and consists of piecewise planar segments The mesh size parameter is the maximal diameter of these segments This maximal diameter is approximately the maximal diameter of the tetrahedra in J that contain the discrete interface i e h hr is approximately the maximal diameter of the tetrahedra in 7 that are close to the interface In Figure 4 1 we illustrate this construction for the two dimensional case Each of the planar segments of I is either a triangle or a quadrilateral The quadrilaterals can formally be divided into two triangles Thus D consists of a set of triangular faces Concerning the maximal distance between I and IT we note the following If we assume II amp n z x ch for all x in a neighbourhood of T which is reasonable for a smooth and piecew
70. ch has its zero level set precisely at the interface T in combination with the Heaviside function H R R Hity fo 05 H 1 fo gt 0 For ease one can take H 0 5 We define p p p2 pi H 9 ud pa u2 um H 6 1 3 Combination of the CSF approach with the level set method leads to the following model for the two phase problem in Q x 0 7 u V u Vp p g div p D u rK rnr divu 0 1 4 Pt Tu Vo 0 together with suitable initial and boundary conditions for u and amp cf section 1 2 3 NS2 combined with transport of a dissolved species NS2 T We consider a two phase flow problem as described above in NS2 We assume that one or both phases contain a dissolved species that is transported due to convection and diffusion and does not adhere to the interface The concentration of this species is denoted by c z t This flow problem can be modeled by the equations 1 4 for the flow variables and a convection diffusion 1 1 ONE AND TWO PHASE FLOW MODELS IN STRONG FORMULATION 11 equation for the concentration c At the interface we need interface conditions for c The first interface condition comes from mass conservation which implies flux continuity The second condition results from a constitutive equation known as Henry s law derived from continuity of chemical potentials at the interface which states that the quotient of the concentration values at both sides of the i
71. contained in U For a sufficiently smooth function g U R the tangential derivative along I is defined by projecting the derivative on the tangent space of I i e Vrg Vg Vg np ny 4 12 42 SPATIAL DISCRETIZATION 35 The Laplace Beltrami operator on I is defined by Arg Vp Vrg It can be shown that Vrg and Arg depend only on values of g on I For vector valued functions f 9 T R we define 3 Arf Arfi Arf Ar fs Vrf Vrg X Vrfi Vrgi i 1 We recall the following basic result from differential geometry Theorem 1 Letidr I R be the identity onT and K k ko the sum of the principal curvatures For all sufficiently smooth vector functions v on T the following holds x vds Ar idr vds J Vr idr Vrv ds 4 13 T T T Let I be a polyhedral approximation of I as described in section 4 2 2 In view of the result in this theorem an obvious choice for fr vr in 4 3 that is used in e g 19 45 46 2 57 is the following Tos vn rf Vr idr Vp Vh ds vg Vp 4 14 D Here idr denotes the identity D R3 i e the coordinate vector on D In 2 it is shown that for piecewise quadratic functions v the result fr vn can easily be determined exactly i e without any further approximation errors Analysis and numerical experiments in 3 yield that for this choice we have Es f r va fr vn UE 4 15 va EV valli and that this bound is
72. ction 4 2 5 First results of this finite element method applied to the Laplace Beltrami equation on a given surface I are presented in 9 The method will also be used for the spatial discretization of the convection diffusion problem 6 1 Clearly for the numerical treatment of this problem we also need a time discretization method Remainder forthcoming 58CHAPTER 6 TWO PHASE FLOW WITH TRANSPORT OF A SURFACTANT NS2 S Chapter 7 Two phase flow with transport of both a dissolved species and a surfactant at the interface Forthcoming 60CHAPTER 7 TWO PHASE FLOW WITH TRANSPORT OF BOTH A DISSOLVED SPECIES AND A Part II Implementation in DROPS 63 In this second part we discuss some implementation issues related to the numerical treated in part I In chapter 8 we describe some fundamental concepts and the most important classes of DROPS In chapter 9 we give a brief introduction to the parallel version 64 Chapter 8 Fundamental concepts and data structures In this chapter important data structures and algorithms implemented in DROPS are presented Figure 8 1 gives an overview of the main components of the software The different modules are arranged in a diagram such that has two levels of structuring namely in vertical and horizontal direction o E fe o o N o lectors E EP oarse Matrices s iz o ON EN D j E Ensight Matl
73. d in TransportP1CL class file poisson transport2phase h e P1 finite element space for the concentration c and the transformed concentration cf chapter 5 The matrices and vectors are assembled in the member functions SetupInstat System e the 6 scheme for the time integration cf section 5 3 in the function DoStep e the transformation from c to and vice versa with the functions c2ct and ct2c 111 112 CHAPTER 13 A TWO PHASE FLOW PROBLEM WITH MASS TRANSPORT 13 2 Implementation The program is almost identical to that in the previous chapter with the following additional components for the mass transport part e The header file for TransportP1CL class include poisson transport2phase h e The functions Initialcneg Initialcpos to describe the initial values of the concentra tion in and Q2 respectively These values are read from the parameter file in the main function e The boundary conditions for the concentration typedef DROPS BndDataCL lt gt cBndDataCL typedef cBndDataCL bnd val fun c_bnd_val_fun const DROPS BndCondT c_bc 6 DROPS OutflowBC DROPS OutflowBC DROPS OutflowBC DROPS OutflowBC DROPS OutflowBC DROPS OutflowBC const c_bnd_val_fun c_bfun 6 0 0 0 0 0 0 e In the function main we assign the parameters of the transport problem which are read from the parameter file to the global variables needed to construct the TransportP1CL class ini c 0 C cF
74. d concentration we derive a weak formulation Define D 9 Cul D 9 For simplicity we assume homogeneous zero boundary conditions for on Q Furthermore for a given function o we assume an initial condition 0 o The weak formulation is as follows Determine t H4 0 such that for all t 0 T CHOTE v Cat u VEE v D P Vet Vv v v e Hy 5 5 Here c v fog c x v r dx denotes the L scalar product Note that in our application the velocity field u is time dependent i e u u t A solution of this problem satisfies c r 0 in trace sense and thus c Cy satisfies the Henry interface condition in 5 4 Furthermore a solution satisfies D t V n r 0 and thus D Ve n r 0 which is the flux continuity condition in 5 3 If the solution it sufficiently smooth then it satisfies the differential equation 0 ap Cult u Vet DAE Cn amp t D amp t Ae t Cn 9 t in the two subdomains Q and Q2 This can be rewritten as Oc at u Vc D 6 Ac in the two subdomains and for all t 0 T and thus we obtain the convection diffusion equation 5 2 This shows that the problem 5 5 is an appropriate weak formulation of the transport problem Note that to obtain the original physical quantity one has to rescale c Cu o l 5 2 Spatial discretization We describe a finite element discretization of the problem in 5 5 Due
75. ducing weight functions a V Ry for nodes and 8 E R for edges of the graph the computational load a vr of the corresponding tetrahedron T and the amount of communication 3 er for the corresponding face F can be described G T V E o B is called a weighted dual graph Figure 9 3 shows a 2D example for a dual graph For a subset V C V we define a V ever 0 v corresponding to the total load of V For a given partitioning m the set Eeu m T1 To E mT Z m 15 9 2 DISTRIBUTION OF WORK LOAD 83 Figure 9 3 Dual graph for 2D triangulation corresponds to the faces forming the processor boundary where communication takes place The graph partitioning problem is given by the following definition Definition 11 Generalized graph partitioning problem For a constant C gt 1 and a given weighted dual graph V E o 3 find a partitioning m V 1 P such that COStcomm M gt B e min e Ecut m and with V m p The graph partitioning problem belongs to the class of NP hard problems in this sense an optimal partitioning cannot be computed efficiently Nevertheless there are a couple of heuristic approaches with polynomial runtime yielding reasonable results For a survey on this topic we refer to 37 We use the package ParMETIS 72 which realizes a parallel multilevel graph partitioning algorithm described in 59 Based on a partitioning m computed by a graph partitioning
76. duction We consider the two phase flow problem described in chapter 12 but now combined with an additional the mass transport equation Thus the model is as in 1 5 The surface tension coefficient 7 is assumed to be constant independent of the concentration of the surfactant The diffusion coefficients and Henry s constant are D 23210 Dz 5 8 10 9 13 1 Cy 0 75 which are much larger than physically realistic ones We use these large values to have a time scale for the transport process that is comparable the one of the flow dynamics Hence one can better observe the change in the concentration profiles while the droplet rises cf figure 13 1 A homogeneous Neumann boundary condition is imposed for the mass transport problem The initial values of the concentrations inside and outside the droplet are 10 and 0 1 respectively Thus at t 0 the Henry interface condition is not satisfied Since the surface tension coefficient is assumed to be independent of the concentration of the third component there is a coupling only in one direction between the two phase flow problem and the transport equation for the third component Thus per time step we first solve the two phase Navier Stokes plus level set equations and then we use the given velocity and interface to solve the transport equation For the two phase flow problem we use the same methods as in chapter 12 For the transport equation we use the following methods which are implemente
77. e elements is necessary due to the fact that parents and children are linked by pointers Consider for example the situation depicted in Figure 9 1 where a parent T and its child T K T are stored on different processors say 1 and 2 Since pointers from one local memory to another are not allowed in a distributed memory setting we have to use a copy to realize this pointer One could store a copy of T on processor 2 to represent the link between T and T as a pointer on processor q If one does not allow such ghost copies all ancestors and descendants of a tetrahedron must be on the same processor This would cause very coarse data granularity poor load balancing and hence low parallel efficiency For each level k and processor p we introduce a set of master elements Mx p C Gr p and a set of ghost elements G p C G p In the formulation of the conditions below we use the conventions K T if status T NoRef and M 44 p 9 We now formalize the conditions on data distribution as follows 9 1 DATA DISTRIBUTION 79 Definition 7 Admissible hierarchical decomposition The distributed hierarchical decom position H is called an admissible hierarchical decomposition if for all k 1 J the following conditions are fulfilled Al Partitioning of G p The sets of masters and ghosts form a disjoint partitioning of Gx p Vp Mx p UGx p Gk p and Mi p Y Gx p 0 A2 Existence Every element from Cj is represe
78. e integra tion a standard 6 scheme which includes Euler backward and the Crank Nicolson method is available A simple method is used for the decoupling of u p 9 and c in each time step Model NS2 S In this model we have a convection diffusion equation on the moving interface T cf 1 6 For the spatial discretization we use special finite element spaces that are obtained from suitable restriction of the P finite element space corresponding to the tetrahedral triangu lation Furthermore interface adapted quadrature rules are needed The time integration has to be adapted to the special time derivative n that occurs in 1 6 Model NS24 M S Forthcoming spatial discretization time integration multilevel P P FE 0 scheme tetrahedral fractional step hierarchy scheme local refinement quadrature operator splitting local coarsening I P PX P54 SDFEM for mass conservation CHAPTER 1 INTRODUCTION re initialization of I Th discretization of fr special quadrature l I P4 transf for c 0 scheme for c PX for c l l FE space on I5 discr of time quadrature on T derivative tn l NS2 M S treatment of cou iterative solvers plings linearization u p fully coupled PCG MG inexact Uzawa Picard iteration for Schur compl precond linearization GMRES GCR MINRES SSOR Jacobi l J fixed point for special preconditioners decoupling u p due to jumps
79. e second part we discuss implementation issues such as data structures important classes and routines corresponding to the numerical methods In the third part we present several examples which show how numerical simulations of different one and two phase flow models can be realized in DROPS Three dimensional incompressible two phase flow problems with or without mass transport form the main application area of DROPS The following models fit in the DROPS framework and will be explained in more detail in section 1 1 1 Navier Stokes equations for one phase flow NS1 Our code is used by some compiler manufacturers as a benchmark test for their C compilers e g SUN Microsoft 8 CHAPTER 1 INTRODUCTION 2 Navier Stokes equations for two phase flow NS2 3 NS2 combined with transport of a dissolved species NS2 T 4 NS2 combined with transport of a surfactant at the interface NS2 S 5 NS2 combined with transport of both a dissolved species and a surfactant at the interface NS2 T S In the models 2 5 surface tension forces are taken into account We list some main features of the numerical methods that are used in DROPS for the simulation of these models e A multilevel hierarchy of tetrahedral triangulations is used Local refinement and coarsen ing routines are available e For spatial discretization we apply Finite Element FE techniques based on conforming spaces Special FE spaces suitable for functions tha
80. e that is a polyhedral domain Definition 1 Triangulation A finite collection 7 T of tetrahedra T C is called a triangulation of Q if the following holds 1 vol T gt 0 for all TET 2 Urer T Q 3 int S Nint T 0 for all S T T with S AT Definition 2 Consistency A triangulation 7 is called consistent if the intersection of any two tetrahedra in 7 is either empty a common face a common edge or a common vertex Definition 3 Stability A sequence of triangulations 79 71 75 is called stable if all angles of all tetrahedra in this sequence are uniformly bounded away from zero Definition 4 Refinement For a given tetrahedron T a triangulation K T of T is called a refinement of T if K T gt 2 and any vertex of any tetrahedron T K T is either a vertex or an edge midpoint of T In this case T is called a child of T and T is called the parent of T A triangulation 7 41 is called refinement of a triangulation Tp Z Tk 1 if for every T Tp either T Ty 41 or K T C 7y44 for some refinement K T of T Definition 5 Multilevel triangulation A sequence of consistent triangulations M 76 is called a multilevel triangulation of Q if the following holds 1 For0x k J Tk4 is a refinement of Tk 2 For 0 lt k J TEN Gais T Se Tz 19 Zr 20 CHAPTER 2 HIERARCHY OF TETRAHEDRAL GRIDS Definition 6 Hierarchical decomposition Let M 70 77 be a multilevel trian gula
81. efect correction method with step size control taken from 94 Algorithm 1 fixed point defect correction method with step size control Setw 1 Repeat until desired accuracy 26 CHAPTER 3 NAVIER STOKES EQUATIONS FOR ONE PHASE FLOW NS1 1 Calculate the defects d Ax N x x B y b da i Bx cCc 2 Solve the a discrete Oseen system for the corrections v and q A N x v B q di Bv da with accuracy tol 3 Step size control Calculate the step length parameter 2 2 Kx A N x wkv BT in B 0 yl y q The discrete Oseen system in step 2 can be solved using iterative solvers that are treated in section 3 5 3 with 4 Update x y 3 5 Iterative solvers for large sparse linear systems In this section we discuss iterative solvers available in DROPS that can be used for solving the large sparse linear systems that arise after discretization and linearization of Navier Stokes equations 3 5 1 Iterative solvers for discretized scalar elliptic problems Let Ax b be a large sparse linear system that results from the discretization of a scalar elliptic problem for example a Poisson equation or a convection diffusion equation Basic iterative methods for such a problem are e Jacobi method Gauss Seidel method e Successive overrelaxation method SOR and its symmetric variant SSOR e Conjugate gradient method CG e Krylov subspace methods for nonsymmetric problems GMRES BiCGS
82. el Graph Partitioning and Sparse Matrix Ordering Journal of Parallel and Distributed Computing 48 1 pp 71 95 1998 KASKADE A toolbox for adaptive multilevel codes http www zib de Scisoft kaskade2 R KIMMEL J A SETHIAN Computing geodesic paths on manifolds Proc Natl Acad Sci 95 pp 8431 8435 1998 P KLOUCEK F Rys Stability of the fractional step 0 scheme for the nonstationary Navier Stokes equations SIAM J Numer Anal 31 pp 1312 1335 1994 S LANG Parallele numerische Simulation instationdrer Probleme mit adaptiven Methoden auf unstrukturierten Gittern Ph D thesis University of Stuttgart 2001 S LANG UG A parallel software tool for unstructured adaptive multigrids In Parallel Computing Fundamentals amp Applications Proceedings of the International conference ParCo 99 E H D Hollander J R Joubert F J Peters H Sips eds Imperial College Press Delft 1999 BIBLIOGRAPHY 121 65 Message Passing Interface Forum MPI A Message Passing Interface Standard Interna tional Journal of Supercomputer Applications pp 165 414 1994 66 OpenMP a specification for shared memory parallelization http www openmp org 67 S OsHER J A SETHIAN Fronts propagating with curvature dependent speed algorithms based on Hamilton Jacobi formulations J Comp Phys 79 pp 12 49 1988 68 S OSHER R P FEDKIW Level set methods An overview and some recent results J Com put Phys 169
83. element space for two phase incom pressible flows Accepted for publication in Comp Vis Sci 117 118 14 15 16 17 18 19 20 26 27 28 29 30 BIBLIOGRAPHY C TERBOVEN A SPIEGEL D AN MEY S Gross V REICHELT Experiences with the OpenMP Parallelization of DROPS a Navier Stokes Solver written in C Proceedings of the first international workshop on OpenMP IWOMP 2005 2005 C TERBOVEN A SPIEGEL D AN MEY S GROSS V REICHELT Parallelization of the C Navier Stokes Solver DROPS with OpenMP Proceedings of the international conference on parallel computing ParCo 2005 2005 DROPS wEBPAGE http www igpm rwth aachen de DROPS N ANDERSON A BJ RCK A new high Order Method of Regula Falsi Type for Computing the Root of an Equation BIT 13 pp 253 2064 1973 E BANSCH Numerical methods for the instationary Navier Stokes equations with a free capillary surface Habilitation thesis Albert Ludwigs University Freiburg 1998 E BANSCH Finite element discretization of the Navier Stokes equations with a free capil lary surface Numer Math 88 pp 203 235 2001 R E BANK A H SHERMAN A WEISER Refinement algorithms and data structures for regular local mesh refinement in Scientific computing R Stepleman ed North Holland Amsterdam pp 3 17 1983 R E BANK B D WELFERT H YSERENTANT A class of iterative methods for solving saddle point problems Numer
84. er we typically have a very small time step Due to this the weighting of the mass matrix part in the Oseen problem which scales like 1 At is relatively large compared to the diffusion and convection parts Therefore systems with the matrix A are not very hard to solve Simple preconditioners QA suffice to obtain a reasonably efficient method For the Schur complement preconditioner we can use the same one as for the one phase Navier Stokes equations If one considers a two phase flow problem in which the two phases are gas and liquid then the differences in density and viscosity between the two phases is very large In such a case it might be necessary to modify the Schur complement preconditioner such that it takes these large density and viscosity jumps into account Such a preconditioner was introduced in 7 Chapter 5 Two phase flow with transport of a dissolved species NS2 M 5 1 Weak formulation of the transport equation We recall the strong formulation of the two phase flow problem coupled with a transport equation as described in section 1 1 The two phase flow is modeled by cf 1 4 p 2 u V u Vp p g div u D u rK rnr divu 0 5 1 u V 0 together with suitable initial and boundary conditions The transport equation for the concen tration of a dissolved species is given by cf 1 5 ae u Vc D Ac 5 2 D 6 Ve n 0 at the interface 5 3 C CHC at the interface
85. er Stokes equations for two phase flow NS2 4 1 Weak formulation In this section we present a weak formulation of the two phase flow model 1 4 For simplicity we use homogeneous Dirichlet boundary conditions for u We use the same spaces V Q as for the one phase Navier Stokes problems cf 3 1 We also use V H Q We define the bilinear forms m LQ x LQ SR m u v IEG u vdr a Vx VR a u v fae tr D u D v dx b VxQoR baa adivudz and the trilinear form c VxVxV R c uv w f no a Vv war Q Note that the bilinear forms a and m as well as the trilinear form c depend on A weak formulation of the problem 1 4 is as follows Determine u t V p t Q and t V such that with fr v I rK rnr Vrr r vdr f rKnr Vrr v ds Q I Here r denotes a Dirac distribution and Vr denotes the tangential gradient as defined in 4 12 below Note that the functional fr is dependent and therefore also dependent 3l 32 CHAPTER 4 NAVIER STOKES EQUATIONS FOR TWO PHASE FLOW NS2 4 2 Spatial discretization 4 2 1 Galerkin finite element discretization Let M 79 77 bea multilevel triangulation of With each triangulation 7 0 k lt J we associate a mesh size parameter h hy Let V C V Qn C L2 Q and Vj C V be standard polynomial finite element spaces corresponding to the triangulation Tp For Vj and V we use piecewise quadratics The c
86. es in the menu icon File In the new window Module Browser choose the module drops and check out them Subdirectory drops will be created in your current directory e Go to this new subdirectory drops cd drops 10 2 Compilation of DROPS You may have to change system dependent settings They are contained in the file drops conf Use a text editor for example nedit to change the file drops conf nedit drops conf change the file drops con f Set your architecture and compiler option ARCH to one of the following to one of e ARCH LINUX for gcc compiler e ARCH INTEL for icc Intel C compiler e ARCH SGI for CC compiler in OS Irix e ARCH SUN for CC compiler in OS Solaris and correct the translator call option ARCH CXX and its options option ARCH CXXFLAGS in the file arch ARC H mk conf For example for ARCH LINUX cd arch LINUX nedit mk conf change the file mk con f cd go back into the root directory drops Now that you have changed the files drops conf and arch ARC H mk conf run the make utility at the shell prompt in the current directory with the target dep make dep This will create the corresponding file dependency needed for compiling DROPS 10 3 BUILDING THE BINARY FILE 95 10 3 Building the binary file There are four main subdirectories with the program examples for e a Poisson problem subdirectory POISSON e a one phase Stokes problem subdirectory STOKES
87. evel pidx 0 amp lset and set the indices to their member matrix and vector description objects via the member function SetIdx lset Phi SetIdx lidx Stokes b SetIdx vidx Stokes v SetIdx vidx cplN SetIdx vidx Stokes c SetIdx pidx Stokes p SetIdx pidx Stokes A SetIdx vidx vidx Stokes B SetIdx pidx vidx Stokes prM SetIdx pidx pidx Stokes M SetIdx vidx vidx Stokes N SetIdx vidx vidx The levelset function and velocity are initialized with the functions DistanceFct and Null respectively lset Init DistanceFct Stokes InitVel amp Stokes v Null To export the solution at each time step to the ensight ouput files for visualization we use an object of the class EnsightP2Sol0utCL EnsightP2Sol0utCL ensight MG lidx const string filename C EnsDir C EnsCase const string datgeo filename geo datpr filename pr datvec filename vel datscl filename scl ensight CaseBegin string C EnsCase case c_str C num steps i ensight DescribeGeom Messzelle datgeo true ensight DescribeScalar Levelset datscl true ensight DescribeScalar Pressure datpr true ensight DescribeVector Velocity datvec true The iterative nonlinear and linear and solvers including the corresponding preconditioners for the Navier Stokes problem are similar to those in the one phase flow example in the previous chapter 12 2 IMPLEMENTATION 107 typedef JACPcCL A
88. f SolverAsPreCL lt ASolverT gt APcT APcT Apc Asolver We create the preconditoner Qs with the name BBTpc typedef MinCommPreCL BBTPcT BBTPcT BBTpc 0 B gt Data M gt Data M pr Data The constructor of the class MinCommPreCL requires the pointer A to the upper left matrix A of the block matrix as its first argument This matrix can only be obtained after we have an object of class InstatNavStokesThetaSchemeCL which should be defined later because its template parameters in turn depend on MinCommPreCL Therefore we initialize A with the null pointer and will set the actual pointer to it when the matrix A is available with the function SetupMatrixA We construct the block preconditioner oseenpc for the Oseen problem typedef BlockPreCL lt APcT BBTPcT false OseenPcT OseenPcT oseenpc Apc BBTpc The third template parameter of the class BlockPreCL indicates the type of the preconditioner The value true denotes the diagonal block preconditioner and false stands for the upper triangular one With this preconditoner we define the block preconditoned GCR solver oseensolver typedef GCRSolverCL OseenPcT OseenBaseSolverT OseenBaseSolverT oseensolver0 oseenpc 100 oseen maxiter 1e 4 false kos typedef BlockMatrixSolverCL lt OseenBaseSolverT gt OseenSolverT OseenSolverT oseensolver oseensolverO0 11 3 RESULTS 101 For the nonlinear solver statsolver we use the fixed point defect correction with step size control
89. following should be done e Download the register form fill in and send it to Computing Center RWTH Aachen Uni versity After that you will receive the login account lt User ID gt by TIM and the password e Login by TIM and unlock Hochleistungsrechnen RWTH Aachen e Write a mail to lt terboven rz rwth aachen de gt and ask to record the user lt User ID gt for the DROPS CVS repository e Configure the computer cluster ssh so that not every time the password is demanded Do the following Login at the computer cluster at Computing Center RWTH Aachen University ssh User IDOcluster rz rwth aachen de Create the ssh subdirectory mkdir ssh Create with the vi text editor the authorized keys file vi authorized keys and copy the contents of the ssh id_rsa pub file from you computer into this file e Configure the DROPS CVS repository on your computer export CVS RSH ssh export CVSROOT lt User ID gt cluster rz rwth aachen de home drops CVS_REPOSITORY If the registration is completed successfully perform the following actions 1This address may change In this case ask the maintainers of DROPS about the current contact person 93 94 CHAPTER 10 HOW TO GET STARTED e Create a work directory for installation mkdir your installation directory and checkout drops from DROPS CVS repository by cvs co drops Alternatively you can use a GUI frontend to CVS such as tkevs tkcvs Choose the option Browse Modul
90. funtion we also need a re initialization method A further issue is the discretization of the localized surface tension force in 1 4 For this we use a Laplace Beltrami technique In this type of problems due to surface tension the pressure is discontinuous across the interface For an appropriate treatment of this discontinuity we use a special finite element space P X Due to this discontinuity and discontinuities in density and viscosity we need special quadrature rules After application of a time integration rule we obtain a large nonlinear system of algebraic equations in which u p and are coupled We apply an iterative decoupling strategy to split the coupled problem for u p into two subproblems for u p and 6 respectively If in the model we have very large jumps in density and viscosity across the interface as for example in a liquid gas system then in the iterative solvers we need special preconditioners that are robust with respect to variation in the size of these jumps Model NS24 M For the spatial discretization of the convection diffusion equation for the con centration c we use standard P finite elements Due to the interface condition cy Cyco with Cy 1 in 1 5 the concentration is discontinuous across the interface A special jump removing transformation is used to eliminate this discontinuity An alternative approach is to use a P X finite element space instead of Pi for spatial discretization For the tim
91. h the discrete aproximation of the interface Q1 int T Qo pn Q Q1 Due to the fact that I is piecewise planar these subdomains Q n N T have a relatively simple geometric structure cf Figure 4 2 Figure 4 2 Planar intersections of I and T 7 Quadrature rules for numerical integration over these subdomains are available We outline these numerical integration methods In section 3 2 2 we explained quadrature rules of degree two and five on tetrahedra If TAT we can use these rules Otherwise we use the fact that I is piecewise planar with respect to the triangulation 7 We compute the planar intersection between T and the eight children of T Such a child is denoted by T cf Figure 4 2 Now we iterate over these children and distinguish two cases cf Figure 4 2 40 CHAPTER 4 NAVIER STOKES EQUATIONS FOR TWO PHASE FLOW NS2 e The intersection is a triangle this means that the interface divides the tetrahedron into a tetrahedron and a prism The prism can be subdivided into three tetrahedra Hence we can use the quadrature rules implemented for tetrahedra e The intersection is a quadrilateral the interface divides the tetrahedron into two prisms Both prisms can be subdivided into three tetrahedra Hence we can use the quadrature rules implemented for tetrahedra We also have to compute integrals over the approximate interface l For this we use a quadra ture rule on triangles of degree five
92. h xy X and the 2D coordinate in the local reference frame Note that V may be lie on multiple boundary segments For the example in Figure 8 2 a vertex may be part of up to 3 boundary segments EdgeCL Each edge E is linked to the two vertices Vi Vo which are connected by E If the edge is further refined into two sub edges E E then there is also a link to the midpoint vertex Vm Note that Ey Vi Vm and Ea VmVa If E is located on the boundary then it stores the indices j of the boundary segments with xy xy C Xj Note that an edge can be located on at most 2 boundary segments FaceCL Each face F is linked to its neighboring tetrahedra For a boundary face the index j of the corresponding unique boundary segment X is stored A face F may possess up to 4 neighboring tetrahedra This is the case if F is an inner face connecting two tetrahedra T and 75 which are irregularly refined such that F is not subdivided by the corresponding green refinement rule Then there are two green children Tj K T1 and T K T2 also sharing F as a common face TetraCL Each tetrahedron T is linked to its 4 vertices 6 edges and 4 faces If T gt 0 i e T is not stored in the initial triangulation 79 then T is linked to its parent tetrahedron If T is refined then it is also linked to its children 7T K T T stores the integer values mark T the refinement mark and status T the actual refinement rule These are used in the refi
93. he corresponding TetraCL objects There are different kinds of iterators to access the simplices in the multilevel triangulation The MultiGridCL member functions GetTriangTetraBegin L and GetTriangTetraEnd L return iterators to cycle through all tetrahedra TET of a triangulation Similarly the member functions GetAllTetraBegin L and GetAllTetraEnd L can be used to iterate over all L Te U Gj j 0 where the level T of the iterated tetrahedra T is increasing from 0 to L Similar iterators exist for vertices edges and faces as well The iterators are implemented such that a corresponding for loop can be executed by mul tiple OpenMP threads in parallel 66 This allows for faster computations on shared memory machines As the importance and availability of multicore architectures is growing nowadays and will grow further in the future this is a relevant advantage regarding computational efficiency 8 2 Numerical objects vectors and sparse matrices Vectors In DROPS there are two different type of vectors SVectorCL for short vectors with a handful of entries and VectorCL for vectors with a large number of entries Throughout this chapter we assume that indices always start with the number zero C style numbering SVectorCL lt RowsN gt is a template class with template parameter RowsN for vectors x R owsN with a fixed dimension RowsN It is mostly used for storing coordinates For this pur pose we defined the typedefs Poi
94. he finite element spaces Q that are used for the approximation of the pressure variable in the Galerkin discretization 4 3 Note that in the two phase flow problems that we consider due to surface tension the pressure is smooth in the two subdomains Q but discontinuous across the interface T Moreover the interface T or its approximation D is not aligned to the faces in the tetrahedral triangulation Due to this the finite element space Qn has to be chosen carefully 42 SPATIAL DISCRETIZATION 37 To simplify the presentation we assume that I is known In practice the approach explained below is applied with T replaced by T Let V k kez Z 1 n be the set of all vertices in the triangulation 7 We also assume that x ET for all k For m gt 1 let H 04 U Q2 denote the Sobolev space of functions for which ujo H Q i 1 2 holds We use the notation lull ouo IlullZ o Ill o for u H Q1 U Q5 We introduce the standard linear finite element space W Wp v E C Q vlr e Pi forall T 7 For the approximation of functions in H Q1 U Q2 m gt 1 that are discontinuous across T in trace sense the finite element space W is not suitable In general for u H Qi U Q2 one can not expect a better bound than inf u vllo lt evh ullmanunz veW cf 4 Therefore one should not use the space Qar W in the Galerkin discretization 4 3 of the two phase flow problem T
95. he parent of T are performed on the processor on which the master element of T and its parent are stored The first result in Lemma 1 shows that every T H has at most one ghost copy Moreover due to A5 all leaves T 77 have no ghost copies In this sense the overlap of tetrahedra between the processors is small Distribution of numerical data Let z RY a vector and A RY N a sparse matrix The numbering J 1 N is associated with certain degrees of freedom of the hierarchical decomposition H Based on the distributed hierarchical decomposition H we will define a corresponding distribution of the nu merical data x and A For this purpose we first introduce the notion of a domain decomposition Definition 8 Domain decomposition Let H be a hierarchical decomposition and H its admissible distribution among the processors Due to the conditions A2 and A3 every tetra hedron T H can be assigned a unique processor on which T is stored as a master element In other words we have a well defined function master H 1 P that is given by master T p amp Te Myry p For 0 j J and 1 p P we define Tj p T 6 master T p and Q p U T TET p Then for each 0 lt j lt J the sequence 7 1 7 P is a partition of the triangulation Tj due to A2 A3 and is called the domain decomposition of level j corresponding to the admissible hierarchical decomposition H
96. he preconditioning of A or of the Schur complement one uses an inner Krylov subspace iteration Standard methods like GMRES or BiCGSTAB do not allow such variable preconditioners There are however Krylov subspace methods that can handle such variable preconditioners In DROPS two of such methods are available namely GCR Generalized Conjugate Residuals cf 79 and FGMRES Flexible GMRES cf 79 In our applications we often use GCR with a block preconditioner introduced in 42 of the form Q Br A 0 Ade 3 35 The application of Qj to a vector b is defined by a Jacobi preconditioned BICGSTAB solver with the stopping criterion Ill lt tol roll applied to the linear system Ax b For the Schur complement preconditioner Qz we use the diagonally scaled preconditioner from 42 i e Q BMj BP BMj M B BM j B7 3 36 where M is the diagonal of the matrix M The inverse of B7 Mj B is approximated by using a CG method with the same stopping criterion as of Qj Multigrid method The multigrid method used for the generalized Stokes problem cf section 3 5 2 and 5 can be used without changes for the Oseen problem provided the underlying flow problem is not convection dominated i e a relatively low Reynolds number For problems with strong convec tion this multigrid method has to be modified Such multigrid solvers for convection dominated flows are not available in DROPS yet Chapter 4 Navi
97. his result we conclude that the order of approximation of the modified space WT is the same as that of WT if we take a m In the context of our applications m 1 is a natural choice Therefore the criterion 4 24 with 0 and a 1 is used to decide which basis functions are deleted from WE The constant 6 has to be set by the user The default value is 6 0 1 Implementation issues We comment on two implementation issues The basis of the modified XFEM space contains functions that are discontinuous across I Therefore in the process of assembling the matrices in the finite element discretization integrals over T Op have to be treated carefully This is further discussed in section 4 2 5 The modified XFEM space depends on the position of the interface and thus on the level set function Therefore in an instationary two phase flow problem in each time step the XFEM finite element space can be different This causes additional technical difficulties in the implementation For example one needs suitable interpolation procedures to handle pressure unknowns that are deleted or created due to the change of the XFEM space Furthermore in the linear algebra part one has to deal with the varying dimension of the pressure space 4 2 5 Quadrature At several places integrals over subdomains Q OT i 1 2 occur Here T is a tetrahedron and Q the discrete approximation of the subdomain This approximation is given throug
98. hoice of the pressure finite element space is explained in section 4 2 4 The Galerkin discretization leads to the following variational problem Find u t Va prt Qn and h t V such that for t 0 T m up e t Vn alur t Va clun t un t vn vn Pa t m g va fr vn Vvn Vn 4 3 b un t qn 0 Van Qn 4 4 dn e t vn o un t Vonlt vn o 0 Von V4 4 5 The term fr vn is an approximation of fr v with T being an approximation of the interface I These approximation issues are discussed in sections 4 2 2 and 4 2 3 To make this semi discrete problem well posed we need initial conditions uj 0 and 0 The discretization 4 5 of the linear hyperbolic level set equation is not stable It can be stabilized using a standard streamline diffusion technique This streamline diffusion stabilization applied to the level set equation can be seen as a Petrov Galerkin method with trial space Vj and test functions For each tetrahedron T 7 a stabilization parameter r is chosen The test functions are then defined as O amp x vy x ru amp z t Vv x 2 ET For an analysis of the streamline diffusion method and reasonable choices for the stabilization parameter r we refer to 76 This leads to the following stabilized variant of 4 5 gt bn e t un 0 Von t vn rus t Von 0 Von Vs 4 6 Tet Here r denotes the L scalar p
99. ies a vector ColIndex with N integer entries and a vector Val with N entries of type RealT For a row i the indices from RowBegin i 1 to RowBegin i 1 indicate the range in Val where the values of the non zero entries are stored The column indices of the corresponding values are stored in ColIndex As it is tedious task to compute the sparsity pattern stored in RowBegin and ColIndex we use an intermediate storage format called SparseMatBuilderCL lt RealT gt when setting up a new sparse matrix M The SparseMatBuilderCL first collects and accumulates all entries in a std map based data structure After that a call of the member function Build automatically creates the corresponding SparseMatBaseCL object M and deletes the maps afterwards As maps are often too memory consuming we use them only for initializing M When updating M in subsequent steps the sparsity pattern is reused by default i e access to SparseMatBuilderCL entries directly returns the corresponding SparseMatBaseCL entries in Val If the sparsity pattern should not be reused for example when the extended pressure space 9 changed because the interface has moved all matrix entries should be deleted by a call to the member function clear to force a complete initialization of the matrix 8 3 The connnection between grid and unknowns indices As mentioned before we decided to decouple the geometrical data grid from the numerical data matrices vectors This is advantageous
100. ies u of the phases 0 i 1 2 the surface tension coefficient 7 and the vector of gravitational acceleration g 8 5 Tools for spatial discretization In the discretization procedures Setup of the problem classes several sparse matrices representing the discrete differential operators and vectors for the right hand side have to be constructed This is done by iterating over all tetrahedra T 7 where for a single tetrahedron T contributions to the matrix and vector entries are computed These contributions are integrals over T and the integrands are functions which can be defined locally on T e g basis functions or gradients of basis functions Grid functions For representing the integrands and computing the integrals over T we use LocalP1CL and LocalP2CL objects for linear and quadratic functions respectively and quadrature rules Quad2CL Quad5CL exact for polynomials up to degree 2 or 5 respectively All these classes have a template parameter ValT for the function values and are derived from a common base class GridFunctionCL lt ValT PointsN gt This class stores PointsN values of type ValT which are as sociated to distinct nodes in a tetrahedron described by barycentric coordinates BaryCoordCL cf Section 8 2 For a LocalP1CL object these nodes are the 4 vertices of the tetrahedron for a LocalP2CL object the 6 midpoints of the edges are added For the Quad CL objects the nodes are defined by the quadrature points
101. in Figure 9 3 has 20 nodes and 24 edges whereas the reduced dual graph G 7 in Figure 9 4 has only 8 nodes and 9 edges After computing a load balancing partitioning m V 1 P of the reduced dual graph G T for the data migration we use an migration algorithm described in 52 The migration of the tetrahedra is carried out by means of the DDD package After the migration for the new distributed hierarchical decomposition H the property master T m T p holds In 52 it is shown that for an admissible input hierarchical decomposition the distributed hierarchical decomposition after the migration is again admissible Remark 4 Migration of numerical data If a tetrahedron T is moving from one processor to another also certain vector entries corresponding to the degrees of freedom on T have to be migrated The valid migration of numerical data is a delicate task and will not be discussed in this thesis 9 3 CURRENT STATUS AND OUTLOOK 85 9 3 Current status and outlook The parallel refinement algorithm and load balancing strategy described in 1 served as a starting point for the parallelization of DROPS The work on parallelization is done at the Chair of Scientific Computing RWTH Aachen University Currently the parallel version of DROPS can solve the following problem classes 1 Poisson problems 2 Stokes problems 3 Navier Stokes problems 4 Two phase flow problems In all these cases we can handle ada
102. ion MG with MultiGridCL amp MG NS GetMGO MakeInitialTriangulation MG shell width c level f level 100 CHAPTER 11 NAVIER STOKES EQUATIONS FOR A ONE PHASE FLOW For the spatial discretization we define the variables A B M N b c cp1M cp1N as pointers to the corresponding members of the problem NS After setting the velocity and pressure indices vidx and pidx we create the numbering for the velocity and pressure unknowns with NS CreateNumberingVel MG GetLastLevel vidx NS CreateNumberingPr MG GetLastLevel pidx The members NumUnknowns of vidx and pidx now contain the number of velocity and pressure unknowns We set the dimensions of the matrices and vectors by assigning their row and column indices to vidx and or pidx with the routine SetIdx The initial condition for u when t 0 is set to v by NS InitVel v amp Null Then we assemble the matrices and the right hand side vector with NS SetupInstatSystem A B M NS SetupNonlinear N v cplN 0 0 NS SetupInstatRhs b c cplM O b 0 As we use the same vector b for the right hand side and the couplings with the matrix A the sum of them is contained in b We also assemble the mass matrix NS SetupPrMass amp M pr The preconditoner Q4 named Apc is definen by typedef JACPcCL APcPcT APcPcT Apcpc typedef BiCGStabSolverCL lt APcPcT gt ASolverT BiCGStab based APcT ASolverT Asolver Apcpc 300 1e 2 relative true typede
103. ise quadratic 4 then we have dist P P5 max 2 max o x I n 2 lt chp 4 2 3 Discretization of the curvature localized force term We now discuss the treatment of the surface tension term In the weak formulation this localized force is represented by the linear functional fr as in 4 2 If we assume that the surface tension coefficient T is constant this functional takes the form fr v r Kirmes vez c Kap vas 4 11 Q p with v an element from the space of testfunctions Hd Q in our case The case with a variable surface tension coefficient is discussed in Remark 3 Note that fr is dependent and therefore dependent A difficulty in the discretization of this surface tension force is that due to the curvature one has to deal with second derivatives One possible approach often used in finite difference and finite volume discretizations is to use the representation Lt VO K div nr div rap and discretize the term on the right handside In the finite element setting it is possible to avoid the discretization of second derivatives The approximation of the localized surface tension force is based on a Laplace Beltrami characterization of the curvature For this we have to introduce some elementary notions from differential geometry For ease of notation we choose a fixed t and suppress the time dependence throughout the rest of this section Let U be an open subset in R and T a connected C compact hypersurface
104. ix Ap A1 T Ap with P ASY E Aplp p 1 Remark 3 Computation of distributed stiffness matrix For a stiffness matrix A RN N the local distributed matrix A IR coincides with the stiffness matrix corresponding to the subdomain Q p with triangulation 77 p Thus the local matrices M can be set up in dependently by the different processors p 1 P without any communication Furthermore the parallelization of the Setup routines cf Section 8 4 is a trivial task The conversion of a distributed into an accumulated vector is achieved by summing up the vector entries on processor boundaries which requires communication between adjacent processors Obviously the conversion in the other direction is not unique For computing the matrix vector 82 CHAPTER 9 PARALLELIZATION multiplication j A x we use the accumulated storage X as input and obtain the result jp in a distributed fashion P P x T 5 T a Ai Scd Appl T Ne I App p 1 p 1 Yp Hence the computation of the matrix vector multiplication does not require any communication The scalar product of two vectors z j can be computed efficiently if one of them is stored accumulated for example and the other one distributed yz Then the computation of P P P Chee ig Go p Us d p 1 p 1 p 1 only requires the global summation of P real numbers MPI AllReduce 9 2 Distribution of work load In the simulation of a rising bu
105. l angles Moreover neighboring triangles can have very different areas As is shown in 9 optimal discretization bounds hold if 751550 is shape regular for T n gt o shape regularity is not required e Each quadrilateral in F can be subdivided into two triangles Let 77 be the induced set consisting of only triangles and such that Upe RT I Define WE Yn E Cn vale ePi for all T The space wr is the space of continuous functions that are piecewise linear on the triangles of T a Clearly Vr WI holds There are however situations in which Vr z WI e Let amp lt i lt m be the collection of all vertices of all tetrahedra in wp and the nodal linear finite element basis function corresponding to Then VEI is spanned by the functions dilr 1 lt i m These functions however are not necessarily independent In compu tations we use this generating system d lr 1 lt i m for solving the discrete problem 6 7 Properties that are of interest for the numerical solution of the resulting linear sys tem such as conditioning of the mass and stiffness matrix are analyzed in a forthcoming paper 57 e In the implementation of this method one has to compute integrals of the form ven eds J hes for T Fp T T The domain T is either a triangle or a quadrilateral The first integral can be computed exactly For the second one standard quadrature rules can be applied for example the one discussed in se
106. lel efficiency of linear solvers The most time consuming task in solving two phase flow problems is solving multiple linear systems of equation We want to use a large number of processors for solving such problems in order to be able to store such large data sets and to have relatively short computing times Synchronization points will become a bottleneck for the parallel efficiency Therefore new algorithms and modified Krylov subspace solvers will be developed and implemented in DROPS A good option may be to apply so called s step methods in these solvers Furthermore more parallel preconditioners have to be developed and implemented e Load balancing The current version of the load balancing technique tries to a uniformly distribute the tetrahedra of a given level in most cases of the finest trian gulation and b keep the number of faces on the processor boundary as small as possible Both goals are not optimal In the case of a the number of tetrahedra does not necessarily represent correctly the amount of work each processor has to do For example Dirichlet boundary conditions do not create degrees of freedom on all vertices of tetrahedra at the domain boundary and or the XFEM may introduce extra degrees of freedom in a certain tetrahedron With respect to b we note that just considering faces may not be a correct measure for the amount of communication between processors It is known 9 3 CURRENT STATUS AND OUTLOOK 87 that
107. lem has a matrix vector representation of the form A 8 f asanes Gs 30 CHAPTER 3 NAVIER STOKES EQUATIONS FOR ONE PHASE FLOW NS1 where the parameter is proportional to 1 At cf 3 17 This linear system has a saddle point structure but opposite to the discrete Stokes equation in 3 20 the system matrix K in this case is nonsymmetric For this type of linear systems methods similar to the ones used for the Stokes problem in section 3 5 2 can be used We again distinguish three classes of methods inexact Uzawa method preconditioned Krylov subspace solvers and multigrid methods We briefly address these classes of methods Inexact Uzawa method The inexact Uzawa method has an algorithmic structure as in 3 28 For the Oseen problem the Schur complement is nonsymmetric and therefore we do not apply a CG method in 3 29 but for V we use a preconditioned GMRES or GCR method Also the preconditioners have to be modified For Q4 one can apply a Jacobi or Gauss Seidel iteration a multigrid method or a preconditioned Krylov subspace iteration More delicate is the choice of Qs One possibility is given in 3 36 below Preconditioned Krylov subspace methods One can apply a Krylov subspace method directly to the saddle point system 3 34 and combine this with a block preconditioner K of the form 3 33 It is convenient to allow for a precondi tioner that varies per iteration Such a variable preconditioner arises if for t
108. ling in one direction between the two phase flow problem 5 1 and the transport equation 5 2 5 4 This allows the following obvious decoupling strategy Given values for u t amp t and C c t we first apply a time integration step tn tn to the two phase flow model 5 1 For this we can use the methods treated in section 4 3 and 4 4 This results in approximations for u t 1 and tn41 These can be used in a time integration step c tn c tn41 for example the 6 scheme given in 5 8 Chapter 6 Two phase flow with transport of a surfactant NS2 S We assume that in the two phase flow problem there is a species called tenside or surfactant which adheres to the interface and that the concentration of this surfactant in the two phases is so small that it can be neglected in the model The concentration of this surfactant is denoted by S x t x T We introduce the orthogonal projection P I nn n normal on I Correspondingly for x T we have an orthogonal decomposition u x t Pu z t I P u z t up z t u x t The tangential gradient is defined by Vr PV and divp VI Ar divp Vp The transport of surfactants at the interface is modeled by a convection diffusion equation cf 31 58 btn S divr Sur SKu n DrArS 6 1 The derivative 0 5 stands for the time derivative of S along a normal path For the numerical treatment of such a problem we introduced a new techni
109. m close Now we define the triangulation mg as an object of the BrickBuilderCL std string mesh C meshfile delim xQ size t idx while idx mesh find first of delim std string npos mesh idx std istringstream brick_info mesh brick_info gt gt dx gt gt dy gt gt dz gt gt nx gt gt ny gt gt nz if brick_info dx dz std cerr lt lt error while reading geometry information lt lt mesh lt lt Wn return 1 Y DROPS Point3DCL orig px py pz px 0 dx pylil dy pz 2 dz DROPS BrickBuilderCL builder orig px py pz nx ny nz DROPS MultiGridCL mg builder which will be adaptively refined in a zone of width C ref width up to the level C ref flevel around the interface through an object adap of the class AdapTriangCL The routine Make InitialTriang creates the initial triangulation at the time t 0 and uses the distance function DistanceFct DROPS AdapTriangCL adap mg C ref width 0 C ref_flevel adap MakeInitialTriang DistanceFct We assign a boundary condition to each boundary segment DROPS BndCondT bc 6 DROPS StokesVelBndDataCL bnd val fun bnd fun 6 for DROPS BndIdxT i 0 i num bnd i 1 be i DROPS WallBC bnd_fun i amp Null and define an object prob of the problem class InstatNavierStokes2PhaseP2P1CL typedef ZeroFlowCL CoeffT typedef DROPS InstatNavierStokes2PhaseP2P1CL lt CoeffT gt MyStokesCL MyStokesCL prob
110. merical methods in particular the FE methods for spatial discretization are based on the weak formulation of these partial differential equations as explained in more detail in part I In section 1 2 we address the issue of initial and boundary conditions used in our models The brief summary of numerical methods given in section 1 3 is elaborated in part I 1 1 ONE AND TWO PHASE FLOW MODELS IN STRONG FORMULATION 9 1 1 One and two phase flow models in strong formulation In this section we give the partial differential equations in strong formulation which describe the different one and two phase incompressible flow models that can be solved numerically using DROPS We always assume that the computational domain 2 C R is a bounded polyhedral domain 1 Navier Stokes equations for one phase flow NS1 Let u u x t and p x t be the velocity and pressure For these unknowns we consider the standard incompressible Navier Stokes equations u V u Vp pg div D u in Q divu 0Oin Q 1 1 with a strictly positive density p and viscosisty p and the strain tensor D u Vu Vu The vector g is a known external force gravity Initial and boundary conditions corresponding to these Navier Stokes equations are discussed in section 1 2 Remark 1 The standard case is that p and u are strictly positive constants Using div u 0 we then get Aui div uD u uAu u ui Aus 2 Navier Stokes equations for two pha
111. mmon with FIN we use an orthogonal projection as in the initialization phase Let w w and ws be these common vertices and W either the line or plane through w wa and ws Then d P Pwv v if PyveT dr v Pwr Pwv vll Ww 4 28 min lt j lt m d w v w otherwise The value d Pwv is calculated by linear interpolation of the known values d w 1 j m The approximate distance for v ACT is then determined by d v min dr v T T v with V T A FIN z 0 4 29 After all v ACT have been updated we choose Umin argmin c Ac d v This vertex is removed from ACT and added to FIN All the neighbors of Vmin that are in FAR become an element of ACT and the whole process is repeated After we have determined the approximate distance grid function d v we still have to mark if a vertex is inside or outside the zero level set This is done by using the sign of the given level set function gel d v sign 4 d v 4 30 The calculated values d v uniquely determine a piecewise quadratic function on the trian gulation 7 This function is defined to be the re initialization of pus For this function one can construct an approximate zero level set D 7 as explained in section 4 2 2 The re initialization procedure guarantees that po C UTET T holds and in this sense the change in the discrete interface is small Clearly in general we have py 4 2 7 Mass conservation
112. n finite element discretization Let M 79 77 be a multilevel triangulation of With each triangulation 7 0 lt k lt J we associate a mesh size parameter h hg Let Va C V Qr C L Q and V C V be standard 21 22 CHAPTER 3 NAVIER STOKES EQUATIONS FOR ONE PHASE FLOW NS1 polynomial finite element spaces corresponding to the triangulation 7 We assume the pair Vah Qh to be LBB stable The standard pair used in DROPS is Hood Taylor P P i e piecewise quadratic finite elements for the velocity and piecewise linears for the pressure Using such a pair Vp Qj the Galerkin discretization reads Find u t Va pn t Qh with pp t 1 0 such that for all vj Vp and all qn Qa m un e t Va alur t vn up t un t vn b vn Pn t 8 vn b un t qu 0 for all t 0 T Let i lt icw and yy lt i lt k be the standard nodal bases of the finite element spaces Vn Qn and consider the representations 3 2 N un t So uj t amp j a t uit uv t 3 3 j 1 prt M pius PE m t px 3 4 j l Using this the Galerkin discretization 3 2 can be rewritten as where M e RAIN M amp 6j dr A c RN N Au p VEVE dr Bere By waive dx Q N N u ERY N y n uy V amp amp j de Thus we obtain a system of DAEs for the unknown vector functions u t p t Numerical methods for the sol
113. nd then each of them is subdivided into six tetrahedra A strip with width 4 5 1074 that contains the interface is refined three times more resulting in a mesh size of 0 25 mm near the interface The initial grid is showed in Figure 12 3 We consider the time interval 0 0 75 with the time step At is 5 1074 The nonlinear system for the velocity and pressure is solved with a tolerance 5 1071 and the tolerance parameter for the levelset equation is 10 19 Figure 12 1 Rising droplet Initial state 12 3 RESULTS 109 In Figure 12 2 position and velocities for different times are shown Figure 12 2 Dynamics of the rising droplet at t 0 15 0 3 0 45 0 6 0 75 color coding indicates velocity magnitude 110 CHAPTER 12 NAVIER STOKES EQUATIONS FOR A TWO PHASE FLOW A zoom in that illustrates the velocity field close to the droplet is given in Figure 12 3 Figure 12 3 Velocity vectors at time t 0 6 Figure 12 4 shows the size of the y component of velocity at the droplets barycenter for t from 0 2 to 0 6 These results were used for validation It turns out that this rising velocity is in good agreement with results from measurements that are available in the literature 0 0506 0 0504 0 0502 0 0498 0 0496 0 0494 0 04921 0 049 0 2 0 25 03 0 35 04 0 45 05 0 55 06 0 65 Figure 12 4 Rising velocity of the droplet Chapter 13 A two phase flow problem with mass transport 13 1 Intro
114. nement algorithm cf 1 Furthermore each simplex class contains an UnknownHandleCL object which stores the indices of unknowns belonging to this simplex cf Section 8 3 Multilevel triangulation The class MultiGridCL represents a multilevel triangulation M T 77 cf Definition 5 The data structure is based on the corresponding hierarchical decomposition H 60 67 cf Definition 6 The tetrahedra are stored in J 1 lists each one for the hierarchical surplus C of a different level The vertices edges and faces are stored in a similar manner where the level of such a sub simplex S is defined as S min T T H contains S as sub simplex Furthermore MultiGridCL contains a BoundaryCL object in which all boundary segments are stored The MultiGridCL constructor takes a MGBuilderCL object as input argument which creates the initial triangulation 79 MGBuilderCL serves as an abstract base class from which specific classes can be derived For instance the derived class BrickBuilderCL can be used to generate an initial triangulation of a cuboid shaped domain The member function Refine calls the refinement algorithm described in 1 It requires that the tetrahedra T 7 in the input multilevel triangulation are marked for refinement or 68 CHAPTER 8 FUNDAMENTAL CONCEPTS AND DATA STRUCTURES for coarsening This can be achieved by calling the member functions SetRegRefMark or SetRemoveMark of t
115. nt2DCL Point3DCL and BaryCoordCL which are identical to SVectorCL lt 2 gt SVectorCL lt 3 gt and SVectorCL lt 4 gt respectively The data type VectorCL is used for storing vectors R where N is large and may differ from object to object It is a typedef for VectorBaseCL lt double gt The class VectorBaseCL lt RealT gt is a template class for vectors with entry type RealT and is an ancestor of std valarray lt RealT gt Thus VectorCL derives the benefits of the efficient expression template mechanisms available for arithmetical operations involving valarray objects By setting a debug flag DebugNumericC range checking and other debug features can be enabled which are switched off by default for performance reasons Matrices There are two different types of matrices in DROPS SMatrixCL for small matrices and MatrixCL for large sparse matrices The template class SMatrixCL lt RowsN ColsN gt is used for small matrices M with fixed dimensions Sparse matrices are stored in objects of the type SparseMatBaseCL lt RealT gt where RealT indicates the type of the entries For convenience we introduced a typedef MatrixCL for SparseMatBaseCL lt double gt We use the compressed row storage format CSR which is described in the following For a sparse matrix with m rows and N non zero entries SparseMatBaseCL contains a vector RowBegin JRRowsN xColsN 8 3 THE CONNNECTION BETWEEN GRID AND UNKNOWNS INDICES 69 with m 1 integer entr
116. nted as a master element on level k P Gk U Mp p 1 A3 Uniqueness Every element from G7 is represented by at most one master element on level k V pi p2 Mx p1 1 M pa 0 gt pi po A4 Child parent locality A child master element and its parent as master or ghost are stored on the same processor Vp VTEG VTEK T Te Mesi p gt T Gilp A5 Ghosts are parents Ghost elements always have children Vp VT p K T z 0 A6 Ghost children locality A ghost element and its children are stored on the same processor Vp VTeG p K T C Mk p Remark 1 Consider a consistent initial triangulation 79 Go with a nonoverlapping distri bution of the tetrahedra Go p N Go q 0 for all p Z q In this case all tetrahedra can be stored as masters and there are no ghosts Then the distributed hierarchical decomposition A Go 1 Go P is obviously admissible Two elementary results are given in the following lemma Lemma 1 Let H as in 9 1 be a distributed hierarchical decomposition The following holds 1 If the conditions A3 A5 and A6 are satisfied then for any element from Gp there is at most one corresponding ghost element VT G Vpq T EGk p NG aq gt p 4 2 If the conditions A1 A2 A3 A4 and A6 are satisfied then all children of a parent are stored as master elements on one processor VT ECG dp K T C Mk p 80 CHAPTER 9 PARALLELIZATION A proof
117. nterface equals a given constant Henry s constant which depends on the given phases This model is as follows cf 31 98 Two phase flow model 1 4 combined with Oc Bt Tu Vc D Ac D 9 Vc n 20 at the interface c CHC at the interface The diffusion coefficient is piecewise constant D D Da D1 H Q In the interface condition we use the notation c for cjg restricted to the interface The constant Cy gt 0 is given Henry s constant Using C 1 Cy 1 H the Henry interface condition can also be written as Cc 0 The model has to be combined with suitable initial and boundary conditions cf section 1 2 4 NS2 combined with transport of a surfactant at the interface NS2 S We consider a two phase flow problem as described above in NS2 We assume that there is a species called tenside or surfactant which adheres to the interface and that the concentration of this surfactant in the two phases is so small that it can be neglected in the model The concentration of this surfactant is denoted by S z t r LT We introduce the orthogonal projection P I nn n normal on T Correspondingly for x T we have an orthogonal decomposition u x t Pu z t 4 I P u z t up z t u z t The tangential gradient is defined by Vp PV and divp VE Ar divr Vr The two phase fluid dynamics model with transport of surfactants at the interface is as follows cf 31
118. o improve this poor approximation quality we extend the space W by adding functions that can represent discontinuities across I For the definition of this space we first introduce some further notation To simplify this notation we do not express the dependence on h in our notation for example W instead of Wp The nodal basis in W is denoted by wz pez Let Or be the collection of all tetrahedra that are intersected by I i e Op U T 7 TAT z 0 Let Ri L Q L7 Q i 1 2 be restriction operators vo on Riv m 0 on Q X 9 in L sense We introduce subsets of Z for which the corresponding basis functions have a nonzero intersection with T TH keET ae M2 and supply NT Z 0 TI k eT a EQ and supply AT 4 0 Corresponding spaces are defined by WEF span Riv kE TT i 1 2 We introduce the extended finite element space which is defined as WF RW e RoW Another characterization of W is given by w wew eW Thus the standard linear finite element space W Wp is extended by the spaces WE bel 2 Note that v WF is discontinuous across I supp v C Qr and that v x 0 for all z V 1 In 13 the following approximation property of the XFEM space WT is derived Theorem 2 For integers 0 l m lt 2 the following holds m lu vllaun ch ulm au for all ue H Q4U Q3 4 21 v 38 CHAPTER 4 NAVIER STOKES EQUATIONS FOR TWO PHASE FLOW NS2 This res
119. okes problem we have to consider a coupled system for velocity u pressure p and level set function y cf section 4 4 During the implementation it turned out that the coupling and time discretization should be combined in one class as they are closely connected to each other However the different coupling classes all have a similar structure thus we decided to derive them from a base class CouplLevelsetBaseCL which stores common data members and defines a common abstract interface by means of virtual member functions such as DoStep For the two phase Stokes problem the class CouplLevelsetStokesCL represents a coupled one step 0 scheme For the two phase Navier Stokes problem we implemented several methods for example the following classes e CouplLevelsetNavStokes2PhaseCL generalized 0 scheme as explained in sections 4 3 and 4 4 e CouplLsNsBaenschCL coupled fractional step scheme with operator splitting This method generalizes the method explained in section 3 3 3 to the two phase flow problem cf 2 8 7 ITERATIVE SOLVERS AND PRECONDITIONERS 73 In all these methods result in linear systems of equations that have to be solved All these classes have a template parameter SolverT controlling the type of the iterative solver used in each time step 8 7 Iterative solvers and preconditioners For the implementation of iterative solvers we tried to use a software design that accounts for the nested hierarchy of the solution methods
120. on can be removed by a reparametrization step Finally note that the shifting of the level set function to obtain a better mass conservation introduces a new source of discretization errors 4 3 Time integration For the two phase flow problem the time discretization is based on a generalization of the 6 scheme given in section 3 3 1 for the one phase flow Navier Stokes equations This generalized method is not found in the literature and therefore we describe its derivation in detail The need for a generalization has two causes Firstly opposite to the one phase flow problem the mass matrix M is no longer constant but may vary in time Secondly if in the discretization the XFEM space is used then the matrix B is in general also time dependent due to the dynamics of the interface In the sections below we present two derivations of a generalized 0 schema In the first case we allow M to be time dependent but assume that B does not depend on t In the second case we allow both M and B to be time dependent The resulting schemes are very similar 4 3 1 The generalized 0 scheme for a time independent B We first consider the Navier Stokes part in 4 7 4 8 We use the notation G 6 g frn E t fr 6 0 ACO A N G t HHA Then the Navier Stokes equations can be written as ME BTR GU d amp fe B t 0 4 31 or equivalently d O ME BTEL M G 9 8 1 B t 0 4 32 We assume that
121. onal step 0 scheme is based on a subdivision of each time interval n At n 1 At in three subintervals with endpoints n a At n 1 a At n 4 1 At For given u the approximations u unt1 e yt at these endpoints are defined by urtea u N Zu Fi u t F5 u jure 3 8 yrti a 2 unte e E F ut Fo u tt9 JE 3 9 ZQ utt Eu untl a A Fiat punto fett 3 10 Q We use a popular variant of this scheme cf 74 94 with F Fi Fy OF 1 0 F This scheme is second order accurate and strongly A stable for 1 1 2 a 1 2y2 Q Lan y7 2 l a The Navier Stokes DAE system 3 5 can be rewritten in a equivalent system of ODEs by eliminating the pressure which can be interpreted as a Lagrange multiplier Applying the fractional step 0 scheme to this Navier Stokes problem in ODE form and transforming it back to its original DAE form results in Mat aH zen 0 Ad nN a T N u ejnt E B pr g _ 1 AJAT N d 3 11 Bu 0 ani mim inii i 1 0 A on l a ds Na ger a Bi p tl gto _ plante 4 NaH 3 12 Burtlie 0 woes 0 A art N a t art je Blp _ Lg 1 9 Au e Ne 3 13 Bu 0 3 4 LINEARIZATION METHODS 25 Note that in this method we perform in each time interval nAt n 1 At three successive scheme type substeps The nonlinear problems for the pairs i pr te urt a pr 1 9 u 1 p in these three
122. or k 0 1 gt n 1 Compute the level set vector Pra from the linear system EE OHS Eig n 448 Solve the following equations for uj E P 1 1 on nl onti gt n e 3M T oA Bir 1 eat ON Bist DL it 0B Pin rj ont e 80 8 s fr din ale M 11 8 4 49 B ti 0 50 CHAPTER 4 NAVIER STOKES EQUATIONS FOR TWO PHASE FLOW NS2 The nonlinear system 4 49 is very similar to the discrete Navier Stokes equations 3 17 of the one phase flow problem Note however that the matrices in 4 49 are different due to the fact that the viscosity and density are not constant over the whole domain This nonlinear system has the form Ax N x x B y b 4 50 Bx c and can be solved with the fixed point defect correction method explained in section 3 4 This results in saddle point problems of the form lt A Amann a that have a very similar structure as the ones in section 3 5 3 4 5 Iterative solvers As can be seen from section 4 4 the decoupling and linearization method results in two types of linear problems namely an Oseen i e linearized Navier Stokes problem and a linear discrete hyperbolic problem cf 4 51 and 4 48 The latter can be solved by for example a precondi tioned GMRES or BiCGSTAB method cf section 3 5 1 For the Oseen equations one can use the methods discussed in section 3 5 3 In the two phase flow applications that we consider chapt
123. plement of the matrix K is given by S BA B7 The matrix K has a block factorization 7 k 91 A BT AB I N0 S l fi Solving the problem K y fo lent problem by block forward backward substitution yields the equiva 1 Solve Az fi 3 21 2 Solve Sy Bz fo y LM 1 1 3 22 3 Solve Ax z B y 3 23 In the Uzawa method one applies an iterative solver e g CG to the Schur complement system in step 2 Note that the matrix S in this system is symmetric positive semi definite The A systems that occur in each iteration of this method and in the steps 1 and 3 are solved sufficiently accurate using some fast Poisson solver We consider a simple variant of this method in which the exact solves of the A systems are replaced by inexact ones Let Q4 be a preconditioner of A We use this preconditioner in the steps 1 and 3 and also for the approximation of the Schur complement is step 2 For this we introduce the notation BQj BT 3 24 28 CHAPTER 3 NAVIER STOKES EQUATIONS FOR ONE PHASE FLOW NS1 We use a nonlinear approximate inverse of S denoted by W For each w V w is an approxi mation to the solution z of Sz w We assume that w w z azx lz a for all w 3 25 holds with a given tolerance parameter lt 1 In our implementation V is the preconditioned CG method Let x y be a given approximation to the solution x y Note that using the block factorization of K we get
124. ptive grids that change during the time integration Current status In order to be able to solve these problems with DROPS on a distributed memory computer the following methods are parallelized We use the structure in Figure 8 1 to discuss the parallelized algorithms and data structures and the algorithms and data structures that remain unchanged The refinement algorithm and load balancing strategy have been parallelized This implies that the underlying data structure of multilevel tetrahedral grids is distributed The refinement algorithm and the load balancing algorithm have been extended such that these can handle degrees of freedom on hierarchical tetrahedral grids The initial triangulation 79 can be read from a mesh file generated by the mesh generator GAMBIT 44 and is then distributed among all processors The geometry as well as the finite element solutions can be transfered to Ensight 43 Geomview 47 or to VTK format 80 The setup routines for matrices and right hand sides of the discretized PDEs were easy to parallelize cf section 9 1 T he data structures for finite elements quadrature index clases etc of the sequential version of DROPS are reused in the parallel version with only minor changes In an adaptive setting it is essential to be able to perform interpolation operations on finite element functions between different grids Changing the triangulation or performing load bal ancing operations has to be done ca
125. que cf 9 To explain the main idea of this method we first consider a stationary model problem namely the Laplace Beltrami equation on a given fixed interface I In weak form this problem is as follows For given f L T with fp fds 0 determine u H T with fp uds 0 such that voee ds J fvds for all ve H T 6 2 r r The solution u is unique and satisfies u H T with l ullz2 ry ellfllz2qr and a constant c independent of f For the discretization of this problem one needs an approximation I of T We assume that this approximate manifold is constructed as follows Let 7 n gt o be a family of tetrahedral triangulations of a fixed domain 2 C IR that contains T These triangulations are assumed to be regular consistent and stable Take 7 75 1559 and denote the set of tetrahedra that form Ty by S We assume that T is a closed manifold such that e I can be decomposed as In Urer T 6 3 where for each T there is a corresponding tetrahedron Sr 7 with T Spr 1T and meas2 T gt 0 To avoid technical complications we assume that this S7 is unique i e T does not coincide with a face of a tetrahedron in 75 55 56CHAPTER 6 TWO PHASE FLOW WITH TRANSPORT OF A SURFACTANT NS2 S e Each T from the decomposition in 6 3 is planar i e either a triangle or a quadrilateral Note that the construction described in section 4 2 2 results in approximate interfaces that satisfy these
126. refully in order to keep access to degrees of freedom that are required for the transfer interpolation of FE functions to another triangulation In the current parallel version all algorithms to create and change adaptive grids in time can handle Pj and P5 FE functions cf section 8 3 Within the parallel solvers it is neccessary to transform vectors from the distributed storage format into the accumulated one see Definition 9 Efficient routines for such transformation operations have been implemented A large part of the overall computing time is consumed by iterative solvers for linear systems Parallel versions of CG GMRES GCR QMR and BiCGStab are implemented For precondi tioning these solvers can be combined with a parallel Jacobi or a blocked SSOR method as well as with iterative solvers in the form of preconditioners In order to solve Stokes or Oseen type problems one can use a parallel Krylov subspace method Furthermore a parallel version of the inexact Uzawa algorithm is available cf sec tions 3 5 2 3 5 3 86 CHAPTER 9 PARALLELIZATION To treat the nonlinearity in the discretized the Navier Stokes equations we use a parallel version of the fixed point defect correction method with step size control cf algorithm 1 in section 3 4 The time integration methods methods available in the sequential DROPS version are easy to parallelize Such parallel variants are available Parallelization of the fast marching method F
127. rete approximation of the interface T Ey 0 Discretization of the curvature localized force term Modified finite element space for the pressure XFEM Quadrature sa sa Br Be ek eG Se GR UR RO eR Ret Re initialization method for the level set function M ss c nservation mt dus Ed Ren RECON Boe REDE o 4 3 Eimeintegration 5 2 ed seen bed Lcd pa wes Bex e has Mee d Wo uS S 4 3 1 4 3 2 The generalized 0 scheme for a time independent B The generalized 0 scheme for a time dependent B 3 12 12 15 19 21 21 21 21 22 23 23 24 25 25 26 26 27 29 4 CONTENTS 4 9 8 An implicit Euler type of method 2 2 2 nn nn nn 48 4 3 4 The generalized fractional step 0 Sscheme 2 2 22mm nn 49 4 4 Decoupling and linearization ooa 49 4 5 Iterative solvers 222 less sss 50 5 Two phase flow with transport of a dissolved species NS2 M 51 5 1 Weak formulation of the transport equation a 51 5 2 Spatial discretization 222 ss 52 5 2 1 Standard linear finite element space len 53 5 2 2 Nitsche s method combined with XFEM 53 9 3 Lime integration Soare pai wi ge VERUS d xem ae ae BE Sa nahe 53 5 4 Decoupling and linearization ooo oa a 54 6 Two phase flow with transport of a surfactant NS2 S 55 7 Two phase flow with transport of both a dissolved species and a surfactant at the interfa
128. roduct over the domain T Let i lt j lt n Vj ic lt j lt x and xj 1 lt j lt z be nodal bases of V Qr and Vp respectively These bases induce corresponding representations of the finite element functions in vector form Functions up t Vn pn t Qh and d t V can be represented as N un t X uj D amp j W t ut uv t j 1 rl Im BE mit px j l L td PE hlt PLE j 1 42 SPATIAL DISCRETIZATION 33 For n Vn and uj Vp we introduce the following mass and stiffness matrices M n e RV M n iy plbn amp amp da Alan ERYN Alena 5 n 1 t DENDE de Bent By waive da 9 N u e RF N N uj ij onen VE amp de E u E R E uj j M J xj Ga rur Vxi dx TERYT H u E R H uj gt gt oan Vio bra Yxi dz TETk We also need the following vectors Elen e RY g i f 608 6 da fr On ERY fr dn i fr Ej Using this notation we obtain the following equivalent formulation of the coupled system of ordinary differential equations 4 3 4 4 4 6 Find a t RY p t RX and A t R such that for all t 0 T In addition we have initial conditions for d and 4 2 2 Discrete approximation of the interface T T In this section we explain how a polyhedral approximation I of T is constructed The level set equation for signed distan
129. rth gett g fr En pBlp 44 DECOUPLING AND LINEARIZATION 49 This scheme is similar to the one in remark 7 but now the mass matrices M and E are treated explicitly i e evaluated at t instead of t 1 in the Navier Stokes equations the level set function is treated explicitly and in the level set equation the velocity is treated explicitly Due to this per time step there is a decoupling between the Navier Stokes and level set equation 4 3 4 The generalized fractional step 0 scheme We generalize the method from section 3 3 2 to the two phase flow system 4 7 4 9 Forthcoming 4 4 Decoupling and linearization In the two phase flow case besides the nonlinear system for the velocity ti we also have the nonlinear coupling between the flow variables u p and the level set one A time step in the 0 scheme 4 41 is given by 4M 6A j Q ar Ng ayant 4 gpT gn og iE ME 1 02 But 0 a gt n 1 vac A E 6H d o E E w 3 E gnt l PE The nonlinear coupling between i p and o is decoupled by a fixed point iteration in each iteration which results in a linear system for the level set function and a nonlinear system for the velocity and pressure The fixed point iteration reads e Set ur g 1 0 2 le 2n QV gt h ees TL AT RE 5n 1 antl 5n In i e Initialize Uy and with the value d from the previous time step Iterate f
130. s In each time step the linear system for the new concentration is solved with a tolerance 10 19 The concentration at several points in time are shown in Figure 13 1 114 CHAPTER 13 A TWO PHASE FLOW PROBLEM WITH MASS TRANSPORT Figure 13 1 Dynamics of the mass transport Droplet position and concentration profiles at t 0 15 0 3 0 45 0 6 0 75 color coding indicates concentration value Chapter 14 Appendix 14 1 Parameter file DROPS parameter file for simulation of two phase flow rising droplet time stepping Time NumSteps StepSize ThetaLevelset ThetaStokes flow solver Stokes InnerIter OuterIter InnerTol OuterTol StokesMethod NavStokes Nonlinear Scheme Tol Iter Reduction levelset solver Levelset Tol Iter SD CurvDiff VolCorrection Coupling Tol Iter Stab Y 1300 5e 4 T backward Euler 1 backward Euler 1000 200 1e 12 1e 10 1 no effect 1 time integration O FS op splitting 1 theta scheme 5e 10 20 0 01 1e 10 10000 0 1 1e 10 not used a t m 1 1 till convergence 1 Laplace Beltrami Stabilization 115 116 CHAPTER 14 APPENDIX re initialization of levelset function Reparam Freq 3 0 no reparametrization Method 1 0 1 fast marching without with modification of zero Y adaptive refinement AdaptRef 1 Freq 1 FinestLevel 3 Width 0 45e 3 Y materi
131. s these describe inflow conditions or conditions at walls e g no slip Such Dirichlet conditions are of the form u x t up z t for x IND with a given function up If for example Op corresponds to a fixed wall then a no slip boundary condition is given by u x t 0 for x Np On Ny we prescribe natural boundary conditions which are often used to describe outflow conditions These natural boundary conditions are of the form ong Perno on NN 1 7 with ng the outward pointing normal on IN and pes a given function external pressure For the case pez 0 we thus obtain a homogeneous natural boundary condition For the two phase flow model NS2 in addition we have to consider initial and boundary conditions for the level set funtion The initial condition is x 0 o x in which do is given and should be such that x R do x 0 T 0 Moreover o should be an approximate signed distance function to T 0 To make the linear hyperbolic equation u V 0 well posed one needs boundary conditions on the inflow boundary 005 x 00 u ne 0 The function in the NS2 model 1 4 is introduced for the numerical purpose of capturing the interface and has no physical interpretation There are no natural e g physics based boundary conditions for at the inflow boundary We are only interested in values of close to the interface zero level on and is evolved according to u Vo 0 onl
132. se flow NS2 We consider the situation in which the domain C IR contains two different immiscible incom pressible newtonian phases fluid fluid or fluid gas A model problem is a liquid drop contained in a surrounding fluid The time dependent domains which contain the phases are denoted by Q Q t and 92 Q t with A UN2 Q The interface between the two phases 0 1005 is denoted by T t To model the forces at the interface we make the standard assumption that the surface tension balances the jump of the normal stress on the interface i e we have a free boundary condition cn r rKn with n nr the unit normal at the interface pointing from Q4 in Q2 r the surface tension coefficient material parameter K the curvature of T and the stress tensor i e v pI uD u Due to the assumption that the phases are viscous and immiscible the velocity should be continuous across the interface In combination with the conservation laws of mass and momentum this yields the following standard model p 28 u V u Vp pig div ujD u in M for i 1 2 divu 0 in 9 cn 7TKn u n p 0 1 2 Initial conditions and boundary conditions are discussed in section 1 2 10 CHAPTER 1 INTRODUCTION This model for a two phase incompressible flow problem is often used in the literature The effect of the surface tension can be expressed in terms of a localized force at the interface cf the so called continu
133. sening methods are available too These grid 1 3 OVERVIEW OF NUMERICAL METHODS 13 related methods are also used for all other models For spatial discretization we use the standard Hood Taylor P5 P finite element pair We need quadrature rules to evaluate the integrals that occur in the weak formulation of the problem After spatial discretization we obtain a large system of nonlinear ordinary differential equations coupled with algebraic constraints due to divu 0 i e a DAE system Differential Algebraic Equation For this system we use a numerical time integration rule In DROPS a fractional step 0 scheme and a method based on operator splitting are available Per time step such a time integration rule results in a large nonlinear system of algebraic equations in which velocity u and pressure p are coupled For linearization we apply a standard Picard iteration with steplength optimization After lin earization we have a large sparse linear system of algebraic equations that is of saddle point type Several efficient iterative solvers like for example preconditioned CG PCG and multi grid methods MG are available Model NS2 One very important issue for this class of problems is the discretization of the level set equation For this we use P5 finite elements combined with streamline diffusion sta bilization SDFEM Another topic is the approximation of the zero level of this discretization go of TS 5 Related to the level set
134. ss FileBuilderCL which is an ancestor of MGBuilderCL and is passed to the constructor of MultiGridCL It reads the files written out before by a MGSerializationCL object and recreates the corresponding MultiGridCL object In this way a simulation run that has stopped can be restarted from the last time step where a serialized multilevel triangulation was saved to the file system In a first step the geometrical data is deserialized from the file system using the class FileBuilderCL After that the vectors representing the numerical solutions are restored by means of the class ReadEnsightP2SolCL see the subsequent section Visualization For 3D visualization purposes we mainly use the software package Ensight 43 The class EnsightP2SolOutCL writes out the geometrical information tetrahedra and coordinates of the vertices and the numerical solutions un p y evaluated in all P degrees of freedom using a specific Ensight file format This format can also be read by other visualization packages such as ParaView 70 The class ReadEnsightP2SolCL restores the vectors d p and from the files written out by the class EnsightP2Sol0utCL However this only works properly if the multilevel triangulations at the time of storing and restoring are the same There are interfaces to some other visualization tools as well e GeomMGOutCL GeomSol0utCL for visualization of geometry and numerical solution with Geomview 47 e TecPlotSolOutCL TecPlot
135. ssure variable by restricting to the subspace of discrete divergence free velocities We consider three time integration methods namely the theta scheme the fractional step 0 scheme and a fractional step method based on operator splitting These three methods are treated in the following subsections As default in DROPS for one phase Navier Stokes equations the fractional step 0 scheme is applied 3 3 1 The 6 scheme The 6 scheme applied to the initial value problem 3 6 is is given by urew _ yold Ap TOF A 9 F u Of tnew 1 8 F tora 0 0 1 24 CHAPTER 3 NAVIER STOKES EQUATIONS FOR ONE PHASE FLOW NS1 For 0 0 we have the explicit d scheme which is only first order accurate and not A stable The Crank Nicholson scheme 0 i i is second order accurate and A stable but opposite to the fractional step 0 scheme cf section 3 3 2 not strongly A stable When 0 1 we obtain the implicit Euler scheme which is also first order accurate but strongly A stable Application of this method to the Navier Stokes DAE system results in ME 4 Au n N jg BTB gt n 1 0g 1 6 Au N d u g 3 7 Bi v In each time step a nonlinear system of equations for the unknowns u p has to be solved Iterative methods for solving this system are discussed in section 3 4 3 3 2 fractional step 0 scheme For a given decomposition F Fi F gt and a given parameter a 0 3 the fracti
136. substeps have a similar form 3 3 3 Operator splitting Another approach is introduced in 35 and further analyzed in 62 This method is based on a splitting in the subspace of divergence free functions of the operator F M Au N into F M 0Ad and F d M 1 0 Ad N u Application of this operator splitting method to the problem 3 5 results in the following method cf 35 METTE 0Au T Blp e gto _ 1 6 Au N 3 14 Bu 0 on l a_n a on 1 2n4 1 2n4 1 wusga ee 3 15 g 0Au BTp te dntl gntli e on l T n l M SAL 0AU B p g p Ag Nie a 3 16 Ba 0 An important property of this method is that the nonlinearity and incompressibility condition in the Navier Stokes equations are decoupled The problems in 3 14 3 16 are linear Stokes type of equations and the problem in 3 15 consists of a nonlinear elliptic system for the velocity unknowns 3 4 Linearization methods Using an implicit time stepping scheme we obtain a nonlinear system of algebraic equations in each time step As an example we consider the first step 3 11 in the fractional step scheme note that the other two steps have a similar form The nonlinear system is as follows AM 0A 9N u BTpnto 8 aM 1 6 A N u u 3 17 B 0 This nonlinear system has the form Ax N x x B y b 3 18 Bx c and is solved by a fixed point d
137. t are discontinuous across the interface can be used e level set method for interface capturing is used e A special Laplace Beltrami method for the discretization of surface tension forces is im plemented e A sharp interface approach is applied in the sense that there is no smoothing or regu larization of discontinuities or delta functions at the interface e We use implicit time integration methods in which flow variables surface tension forces and the level set function are strongly coupled e For the solution of large sparse linear systems preconditioned Krylov subspace methods can be applied Several efficient preconditioners are available Inexact Uzawa type solvers for saddle point problems are implemented Multigrid solvers can be used e n MPI based parallel version of DROPS is currently developed For input an interface to the GAMBIT package is implemented and for visualization purposes data can be transfered to Ensight Tecplot Geomview DROPS is written in C In chapter 8 some properties of the implementation are treated Currently a parallel version on DROPS is under development More information on paralleliza tion issues is given in chapter 9 In the remainder of this introduction we describe the models 1 5 more precisely and give an overview of the numerical methods For ease of presentation the partial differential equations used in the models are given in the strong formulation section 1 1 The nu
138. th Fluids 24 pp 671 691 1997 78 Rusten T and Winther R A preconditioned iterative method for saddlepoint problems SIAM J Matrix Anal Appl 13 1992 887 904 79 Y SAAD Iterative Methods for Sparse Linear Systems PWS Publishing Company 1996 80 WILLIAM J SCHROEDER LISA S AVILA WILLIAM HOFFMAN Visualizing with VTK A Tutorial IEEE Computer Graphics and Applications 20 5 2000 81 J A SETHIAN A fast marching level set method for monotonically advancing fronts Proc Natl Acad Sci 93 1591 1996 82 J A SETHIAN Theory algorithms and applications of level set methods for propagating interfaces In Acta Numerica Vol 5 pp 309 395 1996 122 83 84 85 86 87 88 CO So 92 93 94 95 96 97 98 BIBLIOGRAPHY J A SETHIAN Level set methods and fast marching methods Cambridge University Press 1999 SFB 540 Model based Experimental Analysis of Kinetic Phenomena in Fluid Multi phase Reactive Systems http www sfb540 rwth aachen de Silvester D Wathen A Fast iterative solution of stabilised stokes systems Part II using general block preconditioners SIAM J Numer Anal 31 1994 1352 1367 MG A parallel multilevel platform for unstructured grids http rcswww urz tu dresden de jstiller projects mg M SUSSMAN P SMEREKA S OSHER A level set approach for computing solutions to incompressible two phase flow
139. tion of Q For every tetrahedron T M a unique level number T is defined by KT min k T 7 The set Gk C Th Gr Tea LT k is called the hierarchical surplus on level k k 0 J Note that Go 7o Gk Tr Tr for k 1 J The sequence H Oo 07 is called the hierarchical decomposition of M Note that the multilevel triangulation M can be reconstructed from its hierarchical decomposition Remark 2 Let M bea multilevel triangulation and V 0 k J be the corresponding finite element spaces of continuous functions p C Q such that pr P for all T T q gt 1 The refinement property 1 in definition 5 implies nestedness of these finite element spaces Vk C Vi j4 Now assume that based on some error indicator certain tetrahedra in the finest triangulation 7T are marked for refinement In many refinement algorithms one then modifies the finest triangulation 7 resulting in a new one 77 4 Using such a strategy which we call a one level method the new sequence 70 7741 is in general not a multilevel triangulation because the nestedness property 1 in definition 5 does not hold We also note that when using such a method it is difficult to implement a reasonable coarsening strategy In multilevel refinement algorithms the whole sequence M is used and as output one obtains a sequence M 19 77 with Tj To and J J 1 J J 1 In general one has T Z Tp for k gt
140. to the fact that the diffusion coefficient D amp is discontinuous across the interface I the solution has low regularity across I From this and from the fact that the interface is not aligned with the faces in the triangulation is follows that standard finite elements will not have optimal approximation quality for this type of problem cf 55 An alternative is to use the XFEM method explained in section 4 2 4 For the latter method however one needs a modified weak formulation of the transport problem that results in additional technical complications in the implementation In section 5 2 1 we first consider the standard finite element method applied to 5 5 This method is easy to implement but suboptimal In section 5 2 2 it is explained how the XFEM method can be applied to 5 5 resulting in a spatial discretization with a higher order of accuracy than the standard FE method 5 3 TIME INTEGRATION 53 5 2 1 Standard linear finite element space Let V C V Hj Q be the subspace of continuous linear finite elements The Galerkin discretization of 5 5 reads Find amp t V such that for all t 0 7 Cu G t G t vn Culot Tu Ven t vn DEVE Vus V v Vs In addition we need a initial condition 0 that is obtained from interpolation of 0 o The standard nodal basis in V is denoted by i lt j lt K We introduce matrices M ERK Mo CH Wid de H u e R K
141. two phase systems which are of interest for project partners from engineering departments The development of DROPS is mainly carried out at the Chair of Numerical Mathematics RWTH Aachen University Due to the complexity of two phase flow problems we need the ability to perform parallel computations In a tight cooperation the parallelization of DROPS is realized at the Chair of Scientific Computing RWTH Aachen University The DROPS code is written in C Especially the implementation of the iterative solvers heavily uses the object oriented and template programming features of C Some further information including a gallery of simulation examples can be found on the DROPS website 16 More information concerning the numerical methods used and on results of simulations with can be found in a sequence of papers 1 13 In this guide we describe the main features of the DROPS package which has been developed at the Chair of Numerical Mathematics at RWTH Aachen University The guide consists of three parts In this introductory chapter we describe the models of one and two phase flow problems that can be solved numerically using DROPS sections 1 1 and 1 2 and we give an overview of the numerical methods that are available for such simulations In part I we explain these numerical methods Our treatment however is restricted to a rather global description of the methods and for further information references to the literature are given In th
142. ult shows that this space has optimal approximation properties for piecewise smooth functions Therefore the space Qa W is a good choice in the Galerkin discretization of the two phase flow problem This space has been implemented in DROPS We briefly discuss a few specific issues that do not arise if standard finite element spaces like for example W are used Basis used in the XFEM space W Let vx 1 lt k lt n be the standard nodal basis in the finite element space W In 13 it is shown that Pr isksn U Uk kezr U Rave err forms a basis of the extended finite element space WT Thus there is a unique representation v2 V Beet X BY Rive Y BP Rode ve Ww 4 22 k 1 kezr kei This basis is used in the implementation Concerning stability of this basis we note the fol lowing Let Mr be the mass matrix of this basis with respect to the L scalar product and Dr diag Mr In 13 it is proved that the matrix Dr Mr has a spectral condition number that is unformly w r t the mesh size h bounded Moreover the constants that occur in the spectral condition number bounds are also independent of the supports of the basis functions Ripk k Zr In other words a simple scaling is sufficient to control the stability in L of the basis functions with very small supports Modified XFEM space delete functions with very small supports In general there are basis functions Rv WF with very small support in the sense that
143. um surface force CSF model 32 38 This localized force is given by fr TK pnr Here r is a Dirac function with support on I Its action on a smooth test function w is given by f ip Gr x de I s ds This localization approach can be combined with the level set method for capturing the unknown interface We outline the main idea for a detailed treatment we refer to the literature 38 The level set function denoted by x t is a scalar function At the initial time t 0 we assume a function do x such that olx lt 0 for x Q1 0 olx gt 0 for x Q2 0 o x 0 for x I 0 It is desirable to have the level set function as an approximate signed distance function For the evolution of the interface we consider the trace x t of a single particle x 0 N over time A particle on the interface remains on the interface for all time i e for all z 0 I 0 and all t gt 0 we have x t T t This is equivalent to the condition x t t 0 t gt 0 which we extend to the whole domain as a t t x 0 0 for all z 0 Q and all t gt 0 By differentiating this condition with respect to t we obtain V z t z 0 The displacement of a particle coincides with the velocity field i e x u holds Hence one obtains the first order differential equation u V 0 for t gt 0 and EQ The jumps in the coefficients p and u can be described using the level set function whi
144. urPc Stokes solver typedef InexactUzawaCL lt ULPcT SchurPcT gt StokesSolverT StokesSolverT StokesSolver ULPc SchurPc Hence the object StokesSolver represents an inexact Uzawa method For Q4 we chose some iterations of an SSOR preconditioned CG method ULPc applied to the upper left block of the saddle point matrix The Schur complement preconditioner Qs is given by SchurPc We emphasize that there is a conceptual difference between solver objects and preconditioner objects Solver classes are derived from a common base class SolverBaseCL storing the tolerance and the maximum number of iterations i e the stopping criterion as well as the norm of the residual and number of iterations used after the last execution of the solver Each solver class comprises a member function Solve calling the routine of the iterative solver for a given initial guess In contrast each preconditioner class contains the analogon Apply calling the preconditoner for the initial guess 0 74 CHAPTER 8 FUNDAMENTAL CONCEPTS AND DATA STRUCTURES In the following we list the most important solvers and preconditioners available from the DROPS solver toolbox Solvers Navier Stokes linearization methods cf section 3 4 e FixedPtDefectCorrCL Algorithm 1 with step length w 1 e AdaptFixedPtDefectCorrCL Algorithm 1 with step length w as in 3 19 Both are template classes where the template parameter SolverT determines the type
145. uracy for the Oseen solver is increased by a factor 10 after each iteration of the fixed point method The preconditioner for GCR is as in 3 35 This block preconditioner is implemented in the class BlockPreCL in num stokessolver h To apply Qj to a vector b we use a Jacobi preconditioned BICGSTAB solver For Q we use the preconditioner as in 3 36 97 98 CHAPTER 11 NAVIER STOKES EQUATIONS FOR A ONE PHASE FLOW 11 2 Implementation In this section we describe how the problem and the methods outlined in the previous section are implemented in DROPS navstokes drivcav cpp 11 2 1 Input parameters The parameters needed to solve the Navier Stokes problem are the following e parameters for the fixed point iteration fp_tol tolerance fp_maxiter maximal number of iteration deco red the reduction factor which is used to set the new tolerance of the linear solver for the next fixed point iteration e parameters for the time integration theta parameter of the 0 scheme num timestep number of time step time begin and time end begin and end point of the time interval e parameters for the mesh refinement shell width width of the domain to be further refined c level coarsest level f level finest level They will be read from the command line in the main function 11 2 2 Structure of the program We need the following libraries include geom multigrid h include out
146. ution of this system are discussed in section 3 3 3 2 2 Quadrature For the computation of the entries in the matrices in 3 5 integrals have to be determined For this we implemented several quadrature rules In this section we discuss some of these quadrature rules The nodes and weights in these rules are denoted by x and w i 1 m respectively The method Qr is of the form fte dz Qr f vol T Swe i 1 3 3 TIME INTEGRATION 23 We use quadrature rules which have degree of exactness two and five on a tetrahedron These rules are of the form n dx vol T uf f EPn i 1 The nodes are given in barycentric coordinates A quadrature rule of degree two with m 5 is availabe In this method as nodes the four vertices with weights 1 120 are used and the center of the tetrahedron i e coordinates 5 i i i with weight 2 15 Also a method of degree five with m 15 has been implemented The nodes and weights of this method are listed in table 3 1 nodes weights 0 25 0 25 0 25 0 25 16 135 Ai 7 15 34 3 2665 14 15 112900 u 13 3 15 34 Aa 7 15 34 3 2665 14 15 112900 13 3 15 34 As 10 24 15 40 nye Bs 10 2 15 40 Table 3 1 quadrature rule of degree 5 3 3 Time integration Consider an initial value problem of the form d F u f t u 0 uo 3 6 The Navier Stokes system of DAEs in 3 5 takes this form if one elimates the pre
147. velset h The discretization is obtained using a quadratic finite element space with streamline diffusion stabilization The matrices and vectors are assembled using the routines SetupSystem e The Laplace Beltrami method for the discretization of the curvature localized force term cf section 4 2 3 implemented in the function SF ImprovedLaplBeltrami which is called by the member function AccumulateBndIntegral e The discrete approximation of the interface cf section 4 2 2 is implemented in the class InterfacePatchCL 103 104 CHAPTER 12 NAVIER STOKES EQUATIONS FOR A TWO PHASE FLOW e A mass conservation correction implemented in the routine AdjustVolume e The Fast Marching Method for the reparametrization is implemented in the class Fast MarchCL file levelset fastmarch h Its member function Reparam is called by the function ReparamFastMarching of LevelsetP2CL e The time integration is of implicit Euler type cf sections 4 3 3 and is implemented in the class LinThetaScheme2PhaseCL file levelset coupling h e preconditoned GCR method is used to solve the Oseen problem cf section 11 2 12 2 Implementation The rising bubble problem is implemented in levelset risingBubbleAdapt cpp 12 2 1 Input parameters The input parameters are given in a parameter file namely levelset rising bubble butanol water param a text file with the extension param which contains assignments of the form parameter name value
148. w a ghost copy of T is stored on one other processor say q 2 The children of T if they exist are all stored as masters either on processor p or if T has a ghost copy on processor q For T Gk k gt 0 the parent of T or a copy of it is stored on processor p For the multilevel refinement algorithm a crucial point is that for a tetrahedron T one needs information about all children of T cf 27 1 Due to property 2 this information is available on the local processor p or q without communication The first property shows that in a certain sense the overlap of tetrahedra is small In a parallel run of a simulation the computational load has to be distributed more or less equally among the processors Hence an adaptive finite element solver has to be combined with dynamic load balancing and data migration between the processors This is the topic of Section 9 2 The main results concerning the admissible hierarchical decomposition the parallel multilevel refinement method and the load balancing strategy can be summarized as follows e An admissible hierarchical decomposition has the desirable properties 1 small storage overhead and 2 data locality from above This result is given in Section 9 1 e The application of the parallel refinement algorithm to an admissible hierarchical decom position is well defined and results in an admissible hierarchical decomposition This is proved in 1 77 78 CHAPTER 9 PARALLELIZATIO
149. y for a short time interval After this short time a re initialization of is applied cf section 4 2 6 Due to this one can use the level set equation in 1 4 without any boundary conditions on OQ In the model NS2 M in 1 5 one needs in addition initial and boundary conditions for the concentration c The initial condition is c z 0 co x with a given initial concentration co For the boundary conditions the standard ones namely a Dirichlet i e c given on part of 0Q and a Neumann G2 given on part of the boundary condition can be used In model NS2 S in 1 6 one has to prescribe an inititial concentration S x 0 So x x TL for the surfactant If T is a manifold without boundary droplet no boundary conditions for S are needed 1 3 Overview of numerical methods For the numerical simulation of the models 1 5 many numerical methods are implemented in DROPS These methods are explained in part I In this section we present a compact overview of all important methods used in the form of a matrix of methods cf table 1 1 As can be seen from this table we arrange the different methods according to two criteria namely the models for which they are used and the computational method class they belong to To be more specific we briefly address the methods shown in table 1 1 in a row wise order Model NS1 For the one phase Navier Stokes problem we use a multilevel hierarchy of nested tetrahedral meshes Local refinement and coar

DROPS-Package User's guide - Institut für Geometrie und

Contents

Download Pdf Manuals

Related Search

Related Contents

DROPS-Package User's guide - Institut f&uuml;r Geometrie und

Contents

Download Pdf Manuals

Related Search

Related Contents

DROPS-Package User's guide - Institut für Geometrie und