Select Git revision
boundary_and_interface.py
LDDsimulation.py 13.98 KiB
"""
LDD simulation class
"""
import os
import dolfin as df
import numpy as np
import typing as tp
import mshr
import domainPatch as dp
import helpers as hlp
class LDDsimulation(object):
""" Main Class of LDD simulation
LDDsimulation is the main class for the LDD simulation. It contains methods for
setting up and running the LDD simulation for different meshes and domains.
"""
def __init__(self, dimension: int = 2,#
tol: float = None):
hlp.print_once("Init LDD simulation...")
# parameter as set in NSAC_simulation by Lukas Ostrowski
hlp.print_once("Set internal fenics parameter...")
# calculation exactness for the whole simulation.
if tol:
self.tol = tol
else:
self.tol = df.DOLFIN_EPS
# dimension of the simulation. This is by default 2 as I don't intend
# to implement 3D.
self.dim = dimension
#Parameters
# df.parameters["allow_extrapolation"] = True
# df.parameters["refinement_algorithm"] = "plaza_with_parent_facets"
df.parameters["form_compiler"]["quadrature_degree"] = 10
# # To be able to run DG in parallel
# df.parameters["ghost_mode"] = "shared_facet"
# Form compiler options
df.parameters["form_compiler"]["optimize"] = True
df.parameters["form_compiler"]["cpp_optimize"] = True
# df.parameters['std_out_all_processes'] = False
df.set_log_level(df.LogLevel.PROGRESS)
# df.set_log_level(df.LogLevel.INFO)
# CRITICAL = 50, // errors that may lead to data corruption and suchlike
# ERROR = 40, // things that go boom
# WARNING = 30, // things that may go boom later
# INFO = 20, // information of general interest
# PROGRESS = 16, // what's happening (broadly)
# TRACE = 13, // what's happening (in detail)
# DBG = 10 // sundry
hlp.print_once("Set default parameter...")
self.restart_file = None
# # write checkpoints
# self.write_checkpoints = True
# # dont project solution to CG for visualization
# self.project_to_cg = False
# # no sponge zones
# self.sponge = False
## Public variables.
# directory to write files output to.
self.output_dir = None
# list of subdomain meshes. Is populated by self.init_meshes_and_markers()
self.mesh_subdomain = []
# global marker for subdomains and interfaces. Is set by self.init_meshes_and_markers()
self.domain_marker = None
# List of interfaces of class Interface. This is populated by self.init_interfaces
self.interface = []
# The global interface marker. Gets initialised by self.init_interfaces().
self.interface_marker = None
# List of list of indices indicating subdomains adjacent to interfaces.
# self.adjacent_subdomains[i] gives the two indices of subdomains which
# are share the interface self.interface[i]. This gets populated by
# self.init_interfaces()
self.adjacent_subdomains = None
## Private variables
# The number of subdomains are counted by self.init_meshes_and_markers()
self._number_of_subdomains = 0
# variable to check if self.set_parameters() has been called
self._parameters_set = False
# list of objects of class DomainPatch initialised by self._init_subdomains()
self.subdomain = []
def set_parameters(self,#
output_dir: str,#
subdomain_def_points: tp.List[tp.List[df.Point]],#
is_Richards: bool,#
interface_def_points: tp.List[tp.List[df.Point]],#
adjacent_subdomains: tp.List[np.ndarray],#
mesh_resolution: float,#
viscosity: tp.Dict[int, tp.List[float]],#
porosity: tp.Dict[int, float],#
L: tp.Dict[int, tp.List[float]],#
lambda_param: tp.Dict[int, tp.List[float]],#
relative_permeability: tp.Dict[int, tp.Callable[...,None]])-> None:
""" set parameters of an instance of class LDDsimulation"""
self.output_dir = output_dir
self.subdomain_def_points = subdomain_def_points
self.is_Richards = is_Richards
self.mesh_resolution = mesh_resolution
self.interface_def_points = interface_def_points
self.adjacent_subdomains = adjacent_subdomains
self.viscosity = viscosity
self.porosity = porosity
self.L = L
self.lambda_param = lambda_param
self.relative_permeability = relative_permeability
self._parameters_set = True
def initialise(self) -> None:
""" initialise LDD simulation
Public method to call all the init methods of the LDD simulation.
"""
if not self._parameters_set:
raise RuntimeError('Parameters note set!\nYou forgott to run self.set_parameters(**kwds)')
self._init_meshes_and_markers()
self._init_interfaces()
self._init_subdomains()
## Private methods
def _init_meshes_and_markers(self, subdomain_def_points: tp.List[tp.List[df.Point]] = None,#
mesh_resolution: float = None) -> None:
""" Generate meshes and markerfunctions for sudomains.
This method generate meshes and markerfunctions for sudomains and
sets the internal lists of meshes used by the LDDsimulation class.
The following lists and variables are set:
self.mesh_subdomain: List of subdomain meshes. self.mesh_subdomain[0]
is the whole simulation domain, self.mesh_subdomain[i]
is the submesh of subdomain i.
self.domain_marker: Global marker function on self.mesh_subdomain[0]
marking the subdomains by their respective number.
self._number_of_subdomains Internal variable holding the total number of
actual subdomains (global domain not counted).
The method is intended to be used by the self.initialise() method to set
up the simulation but also accepts needed input parameters explicitly for
greater debugging felxibility.
# Parameters
subdomain_def_points #type: tp.List[tp.List[df.Point]] The sub-
domains are passed as list of lists of dolfin
points forming the polygonal boundary of the
subdomains. subdomain_def_points[0] is
to contain the boundary polygon of the whole
domain omega, whereas subdomain_def_points[i]
is supposed to contain the boundary polygone
of subdomain i.
mesh_resolution #type float positive number used by the
mshr mesh generator to influence how coarse
the mesh is. The higher the number, the finer
the mesh.
"""
if not subdomain_def_points:
subdomain_def_points = self.subdomain_def_points
if not mesh_resolution:
mesh_resolution = self.mesh_resolution
### generate global mesh and subdomain meshes. Also count number of subdomains.
# define the global simulation domain polygon first.
domain = mshr.Polygon(subdomain_def_points[0])
for num, subdomain_vertices in enumerate(subdomain_def_points[1:], 1):
subdomain = mshr.Polygon(subdomain_vertices)
domain.set_subdomain(num, subdomain)
# count number of subdomains for further for loops
self._number_of_subdomains = num
# create global mesh. Submeshes still need to be extracted.
mesh = mshr.generate_mesh(domain, mesh_resolution)
self.mesh_subdomain.append(mesh)
# define domainmarker on global mesh and set marker values to domain
# number as assigned in the above for loop.
self.domain_marker = df.MeshFunction('size_t', mesh, mesh.topology().dim(), mesh.domains())
# extract subdomain meshes
for num in range(1, self._number_of_subdomains+1):
self.mesh_subdomain.append(df.SubMesh(mesh, self.domain_marker, num))
def _init_interfaces(self, interface_def_points: tp.List[tp.List[df.Point]] = None,#
adjacent_subdomains: tp.List[np.ndarray] = None) -> None:
""" generate interfaces and interface markerfunctions for sudomains and
set the internal lists used by the LDDsimulation class.
This initialises a list of interfaces and global interface marker functions
used by the LDDsimulation class.
The following lists and variables are set:
self.interface_marker global interface marker
self.interface list containing interface objectso of class
Interface.
The method is intended to be used internally by the self.initialise()
method to set up the simulation but also accepts needed input parameters
explicitly for greater debugging felxibility. In case input parameters
are provided, they override the ones set by self.initialise().
# Parameters
interface_def_points #type: tp.List[tp.List[df.Point]] list of
lists of dolfin points containing the internal
interfaces dividing the domain omega into subdomains.
the indices of this list introduce a numbering
of the interfaces which will determin the list
of neighbours attached to each subdomain.
adjacent_subdomains #type: tp.List[tp.array[int]] List of
list of indices such that adjacent_subdomains[i]
gives the indices of the subdomains that share
the interface defined by the dolfin points in
interface_def_points
"""
if not interface_def_points:
interface_def_points = self.interface_def_points
else:
# overwrite what has been set by init()
self.adjacent_subdomains = adjacent_subdomains
if not adjacent_subdomains:
adjacent_subdomains = self.adjacent_subdomains
else:
# overwrite what has been set by init()
self.adjacent_subdomains = adjacent_subdomains
mesh = self.mesh_subdomain[0]
# define global interface marker. This is a mesh function on edges of the mesh.
self.interface_marker = df.MeshFunction('size_t', mesh, mesh.topology().dim()-1)
self.interface_marker.set_all(0)
for num, vertices in enumerate(interface_def_points):
self.interface.append(dp.Interface(vertices=vertices, #
tol=mesh.hmin()/100,#
internal=True,#
adjacent_subdomains = adjacent_subdomains[num]))#
# the marking number of each interface corresponds to the mathematical
# numbering of the interfaces (starting with 1 not 0 for the first interface).
# The reason is that we want to mark outer boundaries with the value 0.
self.interface[num].mark(self.interface_marker, num+1)
self.interface[num].init_interface_values(self.interface_marker, num+1)
def _init_subdomains(self) -> None:
""" generate list of opjects of class DomainPatch.
"""
is_richards = self.is_Richards
# The list is_richards has the same length as the nuber of subdomains.
# Therefor it is used in the subdomain initialisation loop.
for subdom_num, isR in enumerate(is_richards, 1):
interface_list = self._get_interfaces(subdom_num)
self.subdomain.append(dp.DomainPatch(#
isRichards = isR,#
mesh = self.mesh_subdomain[subdom_num],#
viscosity = self.viscosity[subdom_num],#
porosity = self.porosity[subdom_num],#
interfaces = self.interface,#
has_interface = interface_list,#
L = self.L[subdom_num],#
lambda_param = self.lambda_param[subdom_num],#
relative_permeability = self.relative_permeability[subdom_num]#
))
def _get_interfaces(self, subdomain_index: int, #
adjacent_subdomains: tp.List[np.ndarray] = None) -> np.ndarray:
""" determin all (global) indices of interfaces belonging to subdomain
with index subdomain_index.
This method determins all (global) indices of interfaces belonging to
subdomain with index subdomain_index from adjacent_subdomains.
"""
if not adjacent_subdomains:
adjacent_subdomains = self.adjacent_subdomains
else:
# overwrite what has been set by init()
self.adjacent_subdomains = adjacent_subdomains
interface_of_subdomain =[]
for interface_ind, adj_doms in enumerate(adjacent_subdomains):
# interface [i,j] and [j,i] are not counted seperately in adjacent_subdomains
# therefor we don't care wether subdomain_index == adj_doms[0] or
# subdomain_index == adj_doms[1]. if subdomain_index is one of either,
# the interface belongs to the subdomain with this index.
if np.isin(subdomain_index, adj_doms, assume_unique=True):
interface_of_subdomain.append(interface_ind)
return interface_of_subdomain