diff --git a/LDDsimulation.py b/LDDsimulation.py
index 06aeb1a76bfae38120ee6ac78c44bf096b4105a7..3dfdab81734fb5d5d37136a745615a5cc991500c 100644
--- a/LDDsimulation.py
+++ b/LDDsimulation.py
@@ -90,14 +90,15 @@ class LDDsimulation(object):
def set_parameters(self,#
output_dir: str,#
subdomain_def_points: tp.List[tp.List[df.Point]],#
- is_Richards: bool,#
+ # this is actually an array of bools
+ isRichards: 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]],#
+ porosity: tp.Dict[int, tp.Dict[str, float]],#
+ L: tp.Dict[int, tp.Dict[str, float]],#
+ lambda_param: tp.Dict[int, tp.Dict[str, float]],#
relative_permeability: tp.Dict[int, tp.Dict[str, tp.Callable[...,None]] ],#
saturation: tp.Dict[int, tp.Callable[...,None]],#
timestep: tp.Dict[int, float],#
@@ -106,7 +107,7 @@ class LDDsimulation(object):
""" 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.isRichards = isRichards
self.mesh_resolution = mesh_resolution
self.interface_def_points = interface_def_points
self.adjacent_subdomains = adjacent_subdomains
@@ -242,7 +243,8 @@ class LDDsimulation(object):
self.interface.append(dp.Interface(vertices=vertices, #
tol=mesh.hmin()/100,#
internal=True,#
- adjacent_subdomains = adjacent_subdomains[num]))#
+ adjacent_subdomains = adjacent_subdomains[num],#
+ isRichards = self.isRichards))#
# 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.
@@ -252,13 +254,14 @@ class LDDsimulation(object):
def _init_subdomains(self) -> None:
""" generate list of opjects of class DomainPatch.
"""
- is_richards = self.is_Richards
+ is_richards = self.isRichards
# 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(#
+ subdomain_index = subdom_num,#
isRichards = isR,#
mesh = self.mesh_subdomain[subdom_num],#
viscosity = self.viscosity[subdom_num],#
diff --git a/RR-2-patch-test.py b/RR-2-patch-test.py
index 9cc714754e7500cd7ad5478dd79e58213ff0f92e..41057ba3cf123837deb8842d7336f71c15792454 100755
--- a/RR-2-patch-test.py
+++ b/RR-2-patch-test.py
@@ -42,7 +42,7 @@ interface_def_points = [interface12_vertices]
# adjacent_subdomains[i] contains the indices of the subdomains sharing the
# interface i (i.e. given by interface_def_points[i]).
adjacent_subdomains = [[1,2]]
-is_Richards = [True, True]
+isRichards = [True, True]
Tmax = 2
dt1 = 0.1
@@ -56,8 +56,8 @@ timestep = {
viscosity = {#
# subdom_num : viscosity
- 1 : [1], #
- 2 : [1]
+ 1 : {'wetting' :1}, #
+ 2 : {'wetting' :1}
}
porosity = {#
@@ -68,14 +68,14 @@ porosity = {#
L = {#
# subdom_num : subdomain L for L-scheme
- 1 : [0.25],#
- 2 : [0.25]
+ 1 : {'wetting' :0.25},#
+ 2 : {'wetting' :0.25}
}
lambda_param = {#
# subdom_num : lambda parameter for the L-scheme
- 1 : [4],#
- 2 : [4]
+ 1 : {'wetting' :4},#
+ 2 : {'wetting' :4}
}
## relative permeabilty functions on subdomain 1
@@ -128,7 +128,7 @@ sat_pressure_relationship = {#
simulation = ldd.LDDsimulation()
simulation.set_parameters(output_dir = "",#
subdomain_def_points = subdomain_def_points,#
- is_Richards = is_Richards,#
+ isRichards = isRichards,#
interface_def_points = interface_def_points,#
adjacent_subdomains = adjacent_subdomains,#
mesh_resolution = 2,#
@@ -355,7 +355,7 @@ g2_i = -lambda1*p2_i
for iterations in range(1):
a1 = (L1*u1*v1)*dx1 + (timestep*df.dot(relative_permeability(saturation(p1_i, 1), 1)*df.grad(u1), df.grad(v1)))*dx1 + (timestep*lambda1*u1*v1)*ds1
- rhs1 = (L1*p1_i*v1)*dx1 - ((saturation(p1_i, 1) - saturation(p1_0, 1))*v1)*dx1 + (timestep*(source1 - g1_i)*v1)*ds1
+ rhs1 = (L1*p1_i*v1)*dx1 - ((saturation(p1_i, 1) - saturation(p1_0, 1))*v1)*dx1 +(timestep*source1*v1)*dx1 -(timestep*(g1_i)*v1)*ds1
A1 = df.assemble(a1)
b1 = df.assemble(rhs1)
p_gamma1.apply(A1,b1)
diff --git a/domainPatch.py b/domainPatch.py
index a878315e64992eca9ee31ed6e1b8f19c4cd8debf..b48d556e8d3cf8ba557f5fdca5fd2c64b5bbfd3c 100644
--- a/domainPatch.py
+++ b/domainPatch.py
@@ -19,6 +19,7 @@ import mshr
import numpy as np
# import domainPatch as dp
import typing as tp
+import ufl as fl
# Interface = tp.NewType('Interface', tp.any)
#import domainPatch as domainPatch
@@ -102,15 +103,16 @@ class DomainPatch(df.SubDomain):
### constructor
def __init__(self, #
+ subdomain_index: int,#
isRichards: bool,#
mesh: df.Mesh,#
porosity: float,#
- viscosity: tp.List[float],#
+ viscosity: tp.Dict[str, float],#
# tp.List[Interface] doesn't work, because it hasen't been defined yet.
interfaces: tp.List[object],#
has_interface: tp.List[int],#
- L: tp.List[float],#
- lambda_param: tp.List[float],#
+ L: tp.Dict[str, float],#
+ lambda_param: tp.Dict[str, float],#
relative_permeability: tp.Dict[str, tp.Callable[...,None]],#
saturation: tp.Callable[..., None],#
timestep: float,#
@@ -120,7 +122,8 @@ class DomainPatch(df.SubDomain):
# we need the methods from the parent class, so we call the __init__-method
# from df.SubDomain to access it.
df.SubDomain.__init__(self)
-
+ ### Class Variables/Objects set by the input to the constructor
+ self.subdomain_index = subdomain_index
self.isRichards = isRichards
self.mesh = mesh
self.porosity = porosity
@@ -131,12 +134,21 @@ class DomainPatch(df.SubDomain):
self.lambda_param = lambda_param
self.relative_permeability = relative_permeability
self.saturation = saturation
+ # timestep size, tau in the paper
self.timestep = timestep
- # variable holding the previous iteration.
- self.previous_iteration = None
+ ### Class variables set by methods
# measures are set by self._init_measures()
+ self.iteration_number = 0
self.dx = None
self.ds = None
+ # solution function(s) for the form set by _init_function_space
+ self.pressure = None
+ # variable holding the previous iteration.
+ self.previous_iteration = None
+ # variable holding the pressure of the previous timestep t_n.
+ self.previous_timestep = None
+ # test function(s) for the form set by _init_function_space
+ self.TestFunction = None
self._init_function_space()
self._init_dof_and_vertex_maps()
self._init_boundary_markers()
@@ -145,21 +157,48 @@ class DomainPatch(df.SubDomain):
# END constructor
#### PUBLIC METHODS
- # def governing_form(self, phase: str = 'wetting') -> ufl.object:
- # """ return the governing form of the subdomain"""
- # if self.isRichards:
- # # define measures
- #
- #
- # else:
- # if phase == "wetting":
- # print("governing form wetting phase")
- # elif phase == "nonwetting":
- # print("governing form nonwetting phase")
- # else:
- # raise RuntimeError('missmatch for input parameter phase. \nExpected either phase == wetting or phase == nonwetting ')
- #
- # return form
+ def govering_problem(self, phase: str, iter_num: int = 0) -> tp.Dict[str, fl.Form]:
+ """ return the governing form ant right hand side as a dictionary"""
+ # define measures
+ dx = self.dx
+ ds = self.ds
+ dt = self.timestep
+ if self.isRichards:
+ Lw = self.L['wetting']
+ # this is pw_i in the paper
+ pw = self.pressure['wetting']
+ # this is pw_{i-1} in the paper
+ pw_prev_iter = self.previous_iteration['wetting']
+ pw_prev_timestep = self.previous_timestep['wetting']
+ v = self.testfunction['wetting']
+ Lambda = self.lambda_param['wetting']
+ kw = self.relative_permeability['wetting']
+ S = self.saturation
+ mu_w = self.viscosity['wetting']
+ # we need to have all interfaces in the form
+ interface_forms = []
+ for index in self.has_interface:
+ interface_forms.append((dt*Lambda*pw*v)*ds[index])
+ form = (Lw*pw*v)*dx + (dt/mu_w*df.dot(kw(S(pw_prev_iter))*df.grad(pw), df.grad(v)))*dx + df.sum(interface_forms)
+ # assemble rhs for
+ gi = self._calc_gi(iteration = iter_num)
+ rhs_interface_forms = []
+ for index in self.has_interface:
+ rhs_interface_forms.append((dt*gi*v)*ds[index])
+ rhs = (Lw*pw_prev_iter*v)*dx - ((S(pw_prev_iter) - S(pw_prev_timestep)*v)*dx + dt*source_w*v)*dx - df.sum(rhs_interface_forms)
+ form_and_rhs = {#
+ 'form': form,#
+ 'rhs' : rhs
+ }
+ else:
+ if phase == "wetting":
+ print("governing form wetting phase")
+ elif phase == "nonwetting":
+ print("governing form nonwetting phase")
+ else:
+ raise RuntimeError('missmatch for input parameter phase. \nExpected either phase == wetting or phase == nonwetting ')
+
+ return form
#### PRIVATE METHODS
def _init_function_space(self) -> None:
@@ -168,14 +207,26 @@ class DomainPatch(df.SubDomain):
Note that P1 FEM is hard coded here, as it is assumed in other methods
aswell. This method sets the class variable
self.function_space
+ self.pressure
+ self.TestFunction
"""
if self.isRichards:
self.function_space = {'wetting': df.FunctionSpace(self.mesh, 'P', 1)}
+ self.pressure = df.TrialFunction(self.function_space['wetting'])
+ self.testfunction = df.TestFunction(self.function_space['wetting'])
else:
self.function_space = {#
'wetting' : df.FunctionSpace(self.mesh, 'P', 1),#
'nonwetting' : df.FunctionSpace(self.mesh, 'P', 1)#
}
+ self.pressure = {
+ 'wetting' : df.TrialFunction(self.function_space['wetting']),
+ 'nonwetting' : df.TrialFunction(self.function_space['nonwetting'])
+ }
+ self.testfunction = {
+ 'wetting' : df.TestFunction(self.function_space['wetting']),
+ 'nonwetting' : df.TestFunction(self.function_space['nonwetting'])
+ }
def _init_dof_and_vertex_maps(self) -> None:
""" calculate dof maps and the vertex to parent map.
@@ -440,20 +491,35 @@ class Interface(BoundaryPart):
def __init__(self, #
domain_index: int = None,#
adjacent_subdomains: np.array = None,#
+ # this is an array of bools
+ isRichards: np.array = None,#
**kwds):
# domain_index should be contained in adjacent_subdomains.
- self.domain_index = domain_index
+ # self.domain_index = domain_index
+ self.isRichards = isRichards
self.adjacent_subdomains = adjacent_subdomains
- # _interface_values is a dictionary holding the dof values over the interface
- # for communication. It is initialiesed by the method init_interface_values().
- self._interface_values = None
+ # pressure_values is a dictionary of dictionaries, one for each of the
+ # adjacent subdomains, holding the dof values over the interface of the
+ # pressure (solution) values for communication. It is initialised by the
+ # method init_interface_values().
+ self.pressure_values = None
+ self.pressure_values_prev = None
+ # same as for self.pressure but for the interface decoupling terms gl.
+ self.gl_values = None
+ self.gl_values_prev = None
+ # dictionary holding the current iteration number i for each of the adjacent
+ # subdomains with index in adjacent_subdomains
+ self.current_iteration = dict()
+ self.current_iteration.update({self.adjacent_subdomains[0]: 0})
+ self.current_iteration.update({self.adjacent_subdomains[1]: 0})
+
# make the parent class methods available.
super().__init__(**kwds)
#### Public Methods #####
- def set_domain_index(self, domain_index: int):
- self.domain_index = domain_index
+ # def set_domain_index(self, domain_index: int):
+ # self.domain_index = domain_index
def set_adjacent_subdomains(self, adjacent_subdomains: np.array):
self.adjacent_subdomains = adjacent_subdomains
@@ -500,19 +566,45 @@ class Interface(BoundaryPart):
interface_marker_value: #type int marker value of interface_marker
"""
vertex_indices = self._vertex_indices(interface_marker, interface_marker_value)
- self._interface_values = dict()
- for vert_ind in vertex_indices:
- self._interface_values.update({vert_ind: 0.0})
+ self.pressure_values = dict()
+ #self.pressure_values is a dictionay which contains for each index in
+ # adjacent_subdomains a dictionary for wetting and nonwetting phase with
+ # as a value yet another dictionary holding the interface values. phew.
+ self.pressure_values.update({self.adjacent_subdomains[0]: dict()})
+ self.pressure_values.update({self.adjacent_subdomains[1]: dict()})
+ self.pressure_values_prev = dict()
+ self.pressure_values_prev.update({self.adjacent_subdomains[0]: dict()})
+ self.pressure_values_prev.update({self.adjacent_subdomains[1]: dict()})
+
+ self.gl_values = dict()
+ self.gl_values.update({self.adjacent_subdomains[0]: dict()})
+ self.gl_values.update({self.adjacent_subdomains[1]: dict()})
+ self.gl_values_prev = dict()
+ self.gl_values_prev.update({self.adjacent_subdomains[0]: dict()})
+ self.gl_values_prev.update({self.adjacent_subdomains[1]: dict()})
+ # we initialise all dictionaries to accomodate two phases, to be more
+ # flexible with coupling Richards and Two-phase models.
+ for subdom_ind in self.adjacent_subdomains:
+ for phase in ['wetting', 'nonwetting']:
+ self.pressure_values[subdom_ind].update({phase: dict()})
+ self.pressure_values_prev[subdom_ind].update({phase: dict()})
+ self.gl_values[subdom_ind].update({phase: dict()})
+ self.gl_values_prev[subdom_ind].update({phase: dict()})
+ for vert_ind in vertex_indices:
+ self.pressure_values[subdom_ind][phase].update({vert_ind: 0.0})
+ self.pressure_values_prev[subdom_ind][phase].update({vert_ind: 0.0})
+ self.gl_values[subdom_ind][phase].update({vert_ind: 0.0})
+ self.gl_values_prev[subdom_ind][phase].update({vert_ind: 0.0})
def write_dofs(self, fem_function: df.Function, #
interface_dofs: np.array,#
dof_to_vert_map: np.array,#
local_to_parent_vertex_map: np.array) -> None:
- """ write dofs of fem_function with indices interface_dofs to self._interface_values
+ """ write dofs of fem_function with indices interface_dofs to self.pressure_values
Write the degrees of freedom of the function fem_function with indices
interface_dofs, i.e. the dofs of the interface, to the dictionary of the
- interface self._interface_values.
+ interface self.pressure_values.
# Parameters
fem_function: #type: df.Function, function to save the dofs to
@@ -521,8 +613,8 @@ class Interface(BoundaryPart):
space on the submesh.
local_to_parent_vertex_map: #type: np.array,
"""
- if self._interface_values == None:
- raise RuntimeError("self._interface_values not initiated. You forgot to \
+ if self.pressure_values == None:
+ raise RuntimeError("self.pressure_values not initiated. You forgot to \
run self.init_interface_values")
parent = local_to_parent_vertex_map
@@ -530,18 +622,18 @@ class Interface(BoundaryPart):
for dof_index in interface_dofs:
# fem_function.vector() is orderd by dof and not by vertex numbering of the mesh.
# parent needs a mesh vertex index of the submesh, therefore d2v is needed.
- self._interface_values.update({parent[d2v[dof_index]]: fem_function.vector()[dof_index]})
+ self.pressure_values.update({parent[d2v[dof_index]]: fem_function.vector()[dof_index]})
- print("interface values\n", self._interface_values)
+ print("interface values\n", self.pressure_values)
def read_dofs(self, fem_function: df.Function, #
interface_dofs: np.array,#
dof_to_vert_map: np.array,#
local_to_parent_vertex_map: np.array) -> None:
- """ read dofs off self._interface_values and overwrite the dofs of fem_function
+ """ read dofs off self.pressure_values and overwrite the dofs of fem_function
with indices interface_dofs
- Read the degrees of freedom of the interface self._interface_values and
+ Read the degrees of freedom of the interface self.pressure_values and
overwrite the dofs of the function fem_function with indices
interface_dofs, i.e. update the interface dofs of the function.
.
@@ -555,8 +647,8 @@ class Interface(BoundaryPart):
space on the submesh.
local_to_parent_vertex_map: #type: np.array,
"""
- if self._interface_values == None:
- raise RuntimeError("self._interface_values not initiated. You forgot to \
+ if self.pressure_values == None:
+ raise RuntimeError("self.pressure_values not initiated. You forgot to \
run self.init_interface_values")
parent = local_to_parent_vertex_map
@@ -564,7 +656,7 @@ class Interface(BoundaryPart):
for dof_index in interface_dofs:
# fem_function.vector() is orderd by dof and not by vertex numbering of the mesh.
# parent needs a mesh vertex index of the submesh, therefore d2v is needed.
- fem_function.vector()[dof_index] = self._interface_values[parent[d2v[dof_index]]]
+ fem_function.vector()[dof_index] = self.pressure_values[parent[d2v[dof_index]]]
print("overwritten function \n", fem_function.vector()[:])