-
Notifications
You must be signed in to change notification settings - Fork 14
Guide to developing test problems for proteus mprans
Proteus uses Python modules to specify the input to the air/water model. I will give a brief overview of the organization and philosophy behind the input modules and then give an example.
Proteus is a toolkit for building multiphysics simulators that can support multiple numerical "backends". For this reason the input models are decomposed first into "physics" and "numerics" modules (p- and n-files) and at a higher level these combinations of physics and numerics are combined into a single multiphysics simulation through specification of the coupling in a "split operator" module (so-module). For example the basic incompressible multiphase flow physics and solvers might be specified in two modules
-
twp_navier_stokes_p.py-- the physics, including the coefficients of the equations, domain geometry, boundary, and initial conditions, etc. -
twp_navier_stokes_n.py-- the numerics, including finite element spaces, quadrature rules, solver tolerances, etc.
While these two modules describe the physics and numerics, the simulation is not complete without specifying how the interface between the two phases will be tracked, so additional p/n modules are also included:
-
ls_p.py/ls_n.py-- the level set model for interface tracking -
vof_p.py/vof_n.py-- the volume fraction model for conservation of fluid phase mass (equivalent to volume conservation for incompressible flow) -
redist_p.py/redist_n.py-- the model for re-initializing the level set to a signed distance function -
ls_consrv_p.py/ls_consrv_n.py-- the model for enforcing mass conservation on the level set
Finally, the split operator (so-file) specifies the order in which these models are solved and the manner in which the coupling between the models is controlled (including the "global" time step for all models). For example, for the dambreak problem we use:
-
dambreak_so.py-- the split operator specification for the multiphysics model
Note that the documentation for the input variables that can be specified in the p-,n-, and so-files is at
While not necessary, it is often good practice to define an additional "helper module" that collects information that is used by all models. In the dambreak example this is
-
dambreak.py-- a helper or global module for shared information.
Given this high level structure to the input module collection, we'll provide an example of building up the input for a specific model. We will start with the domain geometry, then proceed to initial and boundary conditions.
To begin with, copy over the input files from something close to what you are interested in and put them under version control:
cd air-water-vv
git fetch origin
git checkout -b my_new_dambreak origin/master
cd 2d
mkdir dambreak2
cd dambreak2
cp ../dambreak/*.py .
git add *.py# "stage" any changes you want to show
git add file1.py file2.py
# "commit" changes to your current branch
# will bring up your default EDITOR
git commit
# changes are now saved locally, share them to GitHub
# use branch name from earlier
git push origin my_new_dambreak
# feel free to repeat this step as you make more changes
# navigate to branch from: https://github.com/erdc-cm/air-water-vv
# should live at: https://github.com/erdc-cm/air-water-vv/tree/my_new_dambreakNext, describe the geometry of your problem. The class hierarchy of Domain objects is described in the proteus.Domain module. For most 2D problems, it is sufficient to use the PlanarStraightLineGraphDomain class, which can represent fairly general domain geometries in a manner appropriate for triangular mesh generators. To make specification of boundary conditions either, we create a mapping of names (strings) to integer tags that partition the boundary segements into sets of segments. The sets of segments must be sufficiently small to decompose the domain into the desired boundary conditions. In the following example, we use a set of tags sufficient for a general 3D box, but we won't use the 'front' and 'back' tags in 2D.
from proteus import Domain
L = (0.584,0.584) #rectangular box extents, used below
#boundary enumeration, useful for referencing boundary sets
boundaries=['left','right','bottom','top','front','back']
boundaryTags=dict([(key,i+1) for (i,key) in enumerate(boundaries)])
#boundaryTags['left'] == 1, etc.
#describe the PSLG
vertices=[[0.0,0.0],#0
[L[0],0.0],#1
[L[0],L[1]],#2
[0.0,L[1]]]#3
vertexFlags=[boundaryTags['bottom'],
boundaryTags['bottom'],
boundaryTags['top'],
boundaryTags['top']]
segments=[[0,1],
[1,2],
[2,3],
[3,0]]
segmentFlags=[boundaryTags['bottom'],
boundaryTags['right'],
boundaryTags['top'],
boundaryTags['left']]
regions=[[0.5*L[0],0.5*L[1]]]
regionFlags=[1]
domain = Domain.PlanarStraightLineGraphDomain(vertices=vertices,
vertexFlags=vertexFlags,
segments=segments,
segmentFlags=segmentFlags,
regions=regions,
regionFlags=regionFlags)For later convenience and verification we can go ahead and generate some text files that can be viewed externally and attached the boundaryTags mapping (dictionary) to the domain object:
domain.boundaryTags = boundaryTags
# the triangle .poly input, can be viewed with triangle's showme app
domain.writePoly("mesh")
# Stanford .ply format, can be viewed with ParaView, etc.
domain.writePLY("mesh")
# Asymptote vector graphics language, requires asymptote application
domain.writeAsymptote("mesh") Since all the equation sets of an air/water model have the same domain, it is convenient to set up this domain in the helper module, domain.py, then, each of the p-files can import this helper file. In other words, put this line near the top of your p-files:
from dambreak import *The basic idea for describing boundary conditions is to write a function that sets the boundary conditions based on the coordinates of boundary points or the boundary tags that we specified above (preferably the latter). In general we can set either the value of the solution variable or the flux of the "conserved quantity". For example, to set the value of the volume fraction along the "top" boundary of the domain to 1.0 (air-satured), you would add to the vof_p.py module code such as
def getDBC_vof(x,flag):
if flag == boundaryTags['top']:
return lambda x,t: 1.0
dirichletConditions = {0:getDBC_vof}The dirichletConditions variable is the important input variable here. It is a dictionary with keys mapping to the solution variables components (only a single entry for vof, three for 2D navier stokes four for 3D navier-stokes, etc.). The function provided to that dictionary could be named anything, but it needs to take as input a point and intenger flax x,flag and return a function of x,t. In this case we use an anonymous function (also known as a lambda).
For flux type boundary conditions the approach is similar, for example if we want no fluid flux into vof equation at all other boundaries the following code would suffice:
def getAFBC_vof(x,flag):
if flag == boundaryTags['top']:# or x[1] >= L[1] - 1.0e-12:
return None
else:
return lambda x,t: 0.0
advectiveFluxBoundaryConditions = {0:getAFBC_vof}
diffusiveFluxBoundaryConditions = {0:{}}Notice that the diffusiveFluxBoundaryConditions dictionary has an empty dictionary inside it. It is empty in this case because the VOF model has no diffusive terms. It the value of the dictionary is itself a dictionary because each equation can have a diffusive term corresponding to each solution component. This is subtle but comes up below when dealing with the Navier-Stokes equations, which have momentum diffusion with respect to each component of velocity.
The recap, the three variables that must be set are the dirichletConditions, advectiveFluxBoundaryConditions, and diffusiveFluxBoundaryconditions dictionaries, the latter being a dictionary-of-dictionaries. The functions that set the boundary conditions, such as getAFBC_vof can be named anything and could be put in the helper module for convenience.
A more complex is example, for the 2D Navier-Stokes boundary conditions is
#constant atmospheric pressure at top
def getDBC_p(x,flag):
if flag == boundaryTags['top']:
return lambda x,t: 0.0
#no lateral slip of air along top (still air)
def getDBC_u(x,flag):
if flag == boundaryTags['top']:
return lambda x,t: 0.0
#allow vertical flow of air along top and slip elsewhere
def getDBC_v(x,flag):
return None
dirichletConditions = {0:getDBC_p,
1:getDBC_u,
2:getDBC_v}
#no flow of mass/volume except at top of tank
def getAFBC_p(x,flag):
if flag != boundaryTags['top']:
return lambda x,t: 0.0
#no flow of u-momentum except at top
def getAFBC_u(x,flag):
if flag != boundaryTags['top']:
return lambda x,t: 0.0
#no flow of v-momentum except at top
def getAFBC_v(x,flag):
if flag != boundaryTags['top']:
return lambda x,t: 0.0
#no viscous stress in u-momentum except at top
def getDFBC_u(x,flag):
if flag != boundaryTags['top']:
return lambda x,t: 0.0
#no viscous stress in v-momentum anywhere (full slip conditions)
def getDFBC_v(x,flag):
return lambda x,t: 0.0
advectiveFluxBoundaryConditions = {0:getAFBC_p,
1:getAFBC_u,
2:getAFBC_v}
diffusiveFluxBoundaryConditions = {0:{},
1:{1:getDFBC_u},
2:{2:getDFBC_v}}Next, we specify the initial conditions for each solution variable. Again, the conditions are specified with a dictionary with keys being the solution components. In the case of initial conditions, the dictionary values are simple objects that define a single member function uOfXT(self,x,t) that defines the solution component for any point in space and optionally time. For example, for a level set function that is initially the signed distance to the surface of a rectangular water column, we could use
#
waterLine_x = 0.146
waterLine_z = 0.292
def signedDistance(x):
phi_x = x[0]-waterLine_x
phi_z = x[1]-waterLine_z
if phi_x < 0.0:
if phi_z < 0.0:
return max(phi_x,phi_z)
else:
return phi_z
else:
if phi_z < 0.0:
return phi_x
else:
return sqrt(phi_x**2 + phi_z**2)
class PerturbedSurface_phi:
def uOfXT(self,x,t):
return signedDistance(x)
initialConditions = {0:PerturbedSurface_phi()}Again, the initialConditions dictionary is the important input variable. The details of implementing an object with a member function uOfXT are up to the programmer. In the actual dambreak example I put the signedDistance function for the water column the helper module because it's useful for several models.
There are many other parameters in the input of a complex model like two-phase flow in general and often also for specific test problems, particularly those corresponding to lab and field experiments. Many of the variable names are self-explanatory or can be understood from context. Going through an example line-by-line and raising questions in the issue tracker for this repository is the best practice for clarifying innput file issues and becoming productive at developing proteus test problems.
Another important resource is the proteus-dev mailing list: