Skip to content
Snippets Groups Projects
Commit 5ca414d7 authored by David Seus's avatar David Seus
Browse files

add new interface data types in Interface.init_interface_values

parent 23d6d7a8
No related branches found
No related tags found
No related merge requests found
...@@ -90,14 +90,15 @@ class LDDsimulation(object): ...@@ -90,14 +90,15 @@ class LDDsimulation(object):
def set_parameters(self,# def set_parameters(self,#
output_dir: str,# output_dir: str,#
subdomain_def_points: tp.List[tp.List[df.Point]],# 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]],# interface_def_points: tp.List[tp.List[df.Point]],#
adjacent_subdomains: tp.List[np.ndarray],# adjacent_subdomains: tp.List[np.ndarray],#
mesh_resolution: float,# mesh_resolution: float,#
viscosity: tp.Dict[int, tp.List[float]],# viscosity: tp.Dict[int, tp.List[float]],#
porosity: tp.Dict[int, float],# porosity: tp.Dict[int, tp.Dict[str, float]],#
L: tp.Dict[int, tp.List[float]],# L: tp.Dict[int, tp.Dict[str, float]],#
lambda_param: tp.Dict[int, tp.List[float]],# lambda_param: tp.Dict[int, tp.Dict[str, float]],#
relative_permeability: tp.Dict[int, tp.Dict[str, tp.Callable[...,None]] ],# relative_permeability: tp.Dict[int, tp.Dict[str, tp.Callable[...,None]] ],#
saturation: tp.Dict[int, tp.Callable[...,None]],# saturation: tp.Dict[int, tp.Callable[...,None]],#
timestep: tp.Dict[int, float],# timestep: tp.Dict[int, float],#
...@@ -106,7 +107,7 @@ class LDDsimulation(object): ...@@ -106,7 +107,7 @@ class LDDsimulation(object):
""" set parameters of an instance of class LDDsimulation""" """ set parameters of an instance of class LDDsimulation"""
self.output_dir = output_dir self.output_dir = output_dir
self.subdomain_def_points = subdomain_def_points self.subdomain_def_points = subdomain_def_points
self.is_Richards = is_Richards self.isRichards = isRichards
self.mesh_resolution = mesh_resolution self.mesh_resolution = mesh_resolution
self.interface_def_points = interface_def_points self.interface_def_points = interface_def_points
self.adjacent_subdomains = adjacent_subdomains self.adjacent_subdomains = adjacent_subdomains
...@@ -242,7 +243,8 @@ class LDDsimulation(object): ...@@ -242,7 +243,8 @@ class LDDsimulation(object):
self.interface.append(dp.Interface(vertices=vertices, # self.interface.append(dp.Interface(vertices=vertices, #
tol=mesh.hmin()/100,# tol=mesh.hmin()/100,#
internal=True,# 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 # the marking number of each interface corresponds to the mathematical
# numbering of the interfaces (starting with 1 not 0 for the first interface). # 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. # The reason is that we want to mark outer boundaries with the value 0.
...@@ -252,13 +254,14 @@ class LDDsimulation(object): ...@@ -252,13 +254,14 @@ class LDDsimulation(object):
def _init_subdomains(self) -> None: def _init_subdomains(self) -> None:
""" generate list of opjects of class DomainPatch. """ 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. # The list is_richards has the same length as the nuber of subdomains.
# Therefor it is used in the subdomain initialisation loop. # Therefor it is used in the subdomain initialisation loop.
for subdom_num, isR in enumerate(is_richards, 1): for subdom_num, isR in enumerate(is_richards, 1):
interface_list = self._get_interfaces(subdom_num) interface_list = self._get_interfaces(subdom_num)
self.subdomain.append(dp.DomainPatch(# self.subdomain.append(dp.DomainPatch(#
subdomain_index = subdom_num,#
isRichards = isR,# isRichards = isR,#
mesh = self.mesh_subdomain[subdom_num],# mesh = self.mesh_subdomain[subdom_num],#
viscosity = self.viscosity[subdom_num],# viscosity = self.viscosity[subdom_num],#
......
...@@ -42,7 +42,7 @@ interface_def_points = [interface12_vertices] ...@@ -42,7 +42,7 @@ interface_def_points = [interface12_vertices]
# adjacent_subdomains[i] contains the indices of the subdomains sharing the # adjacent_subdomains[i] contains the indices of the subdomains sharing the
# interface i (i.e. given by interface_def_points[i]). # interface i (i.e. given by interface_def_points[i]).
adjacent_subdomains = [[1,2]] adjacent_subdomains = [[1,2]]
is_Richards = [True, True] isRichards = [True, True]
Tmax = 2 Tmax = 2
dt1 = 0.1 dt1 = 0.1
...@@ -56,8 +56,8 @@ timestep = { ...@@ -56,8 +56,8 @@ timestep = {
viscosity = {# viscosity = {#
# subdom_num : viscosity # subdom_num : viscosity
1 : [1], # 1 : {'wetting' :1}, #
2 : [1] 2 : {'wetting' :1}
} }
porosity = {# porosity = {#
...@@ -68,14 +68,14 @@ porosity = {# ...@@ -68,14 +68,14 @@ porosity = {#
L = {# L = {#
# subdom_num : subdomain L for L-scheme # subdom_num : subdomain L for L-scheme
1 : [0.25],# 1 : {'wetting' :0.25},#
2 : [0.25] 2 : {'wetting' :0.25}
} }
lambda_param = {# lambda_param = {#
# subdom_num : lambda parameter for the L-scheme # subdom_num : lambda parameter for the L-scheme
1 : [4],# 1 : {'wetting' :4},#
2 : [4] 2 : {'wetting' :4}
} }
## relative permeabilty functions on subdomain 1 ## relative permeabilty functions on subdomain 1
...@@ -128,7 +128,7 @@ sat_pressure_relationship = {# ...@@ -128,7 +128,7 @@ sat_pressure_relationship = {#
simulation = ldd.LDDsimulation() simulation = ldd.LDDsimulation()
simulation.set_parameters(output_dir = "",# simulation.set_parameters(output_dir = "",#
subdomain_def_points = subdomain_def_points,# subdomain_def_points = subdomain_def_points,#
is_Richards = is_Richards,# isRichards = isRichards,#
interface_def_points = interface_def_points,# interface_def_points = interface_def_points,#
adjacent_subdomains = adjacent_subdomains,# adjacent_subdomains = adjacent_subdomains,#
mesh_resolution = 2,# mesh_resolution = 2,#
...@@ -355,7 +355,7 @@ g2_i = -lambda1*p2_i ...@@ -355,7 +355,7 @@ g2_i = -lambda1*p2_i
for iterations in range(1): 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 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) A1 = df.assemble(a1)
b1 = df.assemble(rhs1) b1 = df.assemble(rhs1)
p_gamma1.apply(A1,b1) p_gamma1.apply(A1,b1)
......
...@@ -19,6 +19,7 @@ import mshr ...@@ -19,6 +19,7 @@ import mshr
import numpy as np import numpy as np
# import domainPatch as dp # import domainPatch as dp
import typing as tp import typing as tp
import ufl as fl
# Interface = tp.NewType('Interface', tp.any) # Interface = tp.NewType('Interface', tp.any)
#import domainPatch as domainPatch #import domainPatch as domainPatch
...@@ -102,15 +103,16 @@ class DomainPatch(df.SubDomain): ...@@ -102,15 +103,16 @@ class DomainPatch(df.SubDomain):
### constructor ### constructor
def __init__(self, # def __init__(self, #
subdomain_index: int,#
isRichards: bool,# isRichards: bool,#
mesh: df.Mesh,# mesh: df.Mesh,#
porosity: float,# porosity: float,#
viscosity: tp.List[float],# viscosity: tp.Dict[str, float],#
# tp.List[Interface] doesn't work, because it hasen't been defined yet. # tp.List[Interface] doesn't work, because it hasen't been defined yet.
interfaces: tp.List[object],# interfaces: tp.List[object],#
has_interface: tp.List[int],# has_interface: tp.List[int],#
L: tp.List[float],# L: tp.Dict[str, float],#
lambda_param: tp.List[float],# lambda_param: tp.Dict[str, float],#
relative_permeability: tp.Dict[str, tp.Callable[...,None]],# relative_permeability: tp.Dict[str, tp.Callable[...,None]],#
saturation: tp.Callable[..., None],# saturation: tp.Callable[..., None],#
timestep: float,# timestep: float,#
...@@ -120,7 +122,8 @@ class DomainPatch(df.SubDomain): ...@@ -120,7 +122,8 @@ class DomainPatch(df.SubDomain):
# we need the methods from the parent class, so we call the __init__-method # we need the methods from the parent class, so we call the __init__-method
# from df.SubDomain to access it. # from df.SubDomain to access it.
df.SubDomain.__init__(self) df.SubDomain.__init__(self)
### Class Variables/Objects set by the input to the constructor
self.subdomain_index = subdomain_index
self.isRichards = isRichards self.isRichards = isRichards
self.mesh = mesh self.mesh = mesh
self.porosity = porosity self.porosity = porosity
...@@ -131,12 +134,21 @@ class DomainPatch(df.SubDomain): ...@@ -131,12 +134,21 @@ class DomainPatch(df.SubDomain):
self.lambda_param = lambda_param self.lambda_param = lambda_param
self.relative_permeability = relative_permeability self.relative_permeability = relative_permeability
self.saturation = saturation self.saturation = saturation
# timestep size, tau in the paper
self.timestep = timestep self.timestep = timestep
# variable holding the previous iteration. ### Class variables set by methods
self.previous_iteration = None
# measures are set by self._init_measures() # measures are set by self._init_measures()
self.iteration_number = 0
self.dx = None self.dx = None
self.ds = 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_function_space()
self._init_dof_and_vertex_maps() self._init_dof_and_vertex_maps()
self._init_boundary_markers() self._init_boundary_markers()
...@@ -145,21 +157,48 @@ class DomainPatch(df.SubDomain): ...@@ -145,21 +157,48 @@ class DomainPatch(df.SubDomain):
# END constructor # END constructor
#### PUBLIC METHODS #### PUBLIC METHODS
# def governing_form(self, phase: str = 'wetting') -> ufl.object: def govering_problem(self, phase: str, iter_num: int = 0) -> tp.Dict[str, fl.Form]:
# """ return the governing form of the subdomain""" """ return the governing form ant right hand side as a dictionary"""
# if self.isRichards: # define measures
# # define measures dx = self.dx
# ds = self.ds
# dt = self.timestep
# else: if self.isRichards:
# if phase == "wetting": Lw = self.L['wetting']
# print("governing form wetting phase") # this is pw_i in the paper
# elif phase == "nonwetting": pw = self.pressure['wetting']
# print("governing form nonwetting phase") # this is pw_{i-1} in the paper
# else: pw_prev_iter = self.previous_iteration['wetting']
# raise RuntimeError('missmatch for input parameter phase. \nExpected either phase == wetting or phase == nonwetting ') pw_prev_timestep = self.previous_timestep['wetting']
# v = self.testfunction['wetting']
# return form 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 #### PRIVATE METHODS
def _init_function_space(self) -> None: def _init_function_space(self) -> None:
...@@ -168,14 +207,26 @@ class DomainPatch(df.SubDomain): ...@@ -168,14 +207,26 @@ class DomainPatch(df.SubDomain):
Note that P1 FEM is hard coded here, as it is assumed in other methods Note that P1 FEM is hard coded here, as it is assumed in other methods
aswell. This method sets the class variable aswell. This method sets the class variable
self.function_space self.function_space
self.pressure
self.TestFunction
""" """
if self.isRichards: if self.isRichards:
self.function_space = {'wetting': df.FunctionSpace(self.mesh, 'P', 1)} 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: else:
self.function_space = {# self.function_space = {#
'wetting' : df.FunctionSpace(self.mesh, 'P', 1),# 'wetting' : df.FunctionSpace(self.mesh, 'P', 1),#
'nonwetting' : 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: def _init_dof_and_vertex_maps(self) -> None:
""" calculate dof maps and the vertex to parent map. """ calculate dof maps and the vertex to parent map.
...@@ -440,20 +491,35 @@ class Interface(BoundaryPart): ...@@ -440,20 +491,35 @@ class Interface(BoundaryPart):
def __init__(self, # def __init__(self, #
domain_index: int = None,# domain_index: int = None,#
adjacent_subdomains: np.array = None,# adjacent_subdomains: np.array = None,#
# this is an array of bools
isRichards: np.array = None,#
**kwds): **kwds):
# domain_index should be contained in adjacent_subdomains. # 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 self.adjacent_subdomains = adjacent_subdomains
# _interface_values is a dictionary holding the dof values over the interface # pressure_values is a dictionary of dictionaries, one for each of the
# for communication. It is initialiesed by the method init_interface_values(). # adjacent subdomains, holding the dof values over the interface of the
self._interface_values = None # 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. # make the parent class methods available.
super().__init__(**kwds) super().__init__(**kwds)
#### Public Methods ##### #### Public Methods #####
def set_domain_index(self, domain_index: int): # def set_domain_index(self, domain_index: int):
self.domain_index = domain_index # self.domain_index = domain_index
def set_adjacent_subdomains(self, adjacent_subdomains: np.array): def set_adjacent_subdomains(self, adjacent_subdomains: np.array):
self.adjacent_subdomains = adjacent_subdomains self.adjacent_subdomains = adjacent_subdomains
...@@ -500,19 +566,45 @@ class Interface(BoundaryPart): ...@@ -500,19 +566,45 @@ class Interface(BoundaryPart):
interface_marker_value: #type int marker value of interface_marker interface_marker_value: #type int marker value of interface_marker
""" """
vertex_indices = self._vertex_indices(interface_marker, interface_marker_value) vertex_indices = self._vertex_indices(interface_marker, interface_marker_value)
self._interface_values = dict() self.pressure_values = dict()
for vert_ind in vertex_indices: #self.pressure_values is a dictionay which contains for each index in
self._interface_values.update({vert_ind: 0.0}) # 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, # def write_dofs(self, fem_function: df.Function, #
interface_dofs: np.array,# interface_dofs: np.array,#
dof_to_vert_map: np.array,# dof_to_vert_map: np.array,#
local_to_parent_vertex_map: np.array) -> None: 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 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_dofs, i.e. the dofs of the interface, to the dictionary of the
interface self._interface_values. interface self.pressure_values.
# Parameters # Parameters
fem_function: #type: df.Function, function to save the dofs to fem_function: #type: df.Function, function to save the dofs to
...@@ -521,8 +613,8 @@ class Interface(BoundaryPart): ...@@ -521,8 +613,8 @@ class Interface(BoundaryPart):
space on the submesh. space on the submesh.
local_to_parent_vertex_map: #type: np.array, local_to_parent_vertex_map: #type: np.array,
""" """
if self._interface_values == None: if self.pressure_values == None:
raise RuntimeError("self._interface_values not initiated. You forgot to \ raise RuntimeError("self.pressure_values not initiated. You forgot to \
run self.init_interface_values") run self.init_interface_values")
parent = local_to_parent_vertex_map parent = local_to_parent_vertex_map
...@@ -530,18 +622,18 @@ class Interface(BoundaryPart): ...@@ -530,18 +622,18 @@ class Interface(BoundaryPart):
for dof_index in interface_dofs: for dof_index in interface_dofs:
# fem_function.vector() is orderd by dof and not by vertex numbering of the mesh. # 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. # 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, # def read_dofs(self, fem_function: df.Function, #
interface_dofs: np.array,# interface_dofs: np.array,#
dof_to_vert_map: np.array,# dof_to_vert_map: np.array,#
local_to_parent_vertex_map: np.array) -> None: 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 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 overwrite the dofs of the function fem_function with indices
interface_dofs, i.e. update the interface dofs of the function. interface_dofs, i.e. update the interface dofs of the function.
. .
...@@ -555,8 +647,8 @@ class Interface(BoundaryPart): ...@@ -555,8 +647,8 @@ class Interface(BoundaryPart):
space on the submesh. space on the submesh.
local_to_parent_vertex_map: #type: np.array, local_to_parent_vertex_map: #type: np.array,
""" """
if self._interface_values == None: if self.pressure_values == None:
raise RuntimeError("self._interface_values not initiated. You forgot to \ raise RuntimeError("self.pressure_values not initiated. You forgot to \
run self.init_interface_values") run self.init_interface_values")
parent = local_to_parent_vertex_map parent = local_to_parent_vertex_map
...@@ -564,7 +656,7 @@ class Interface(BoundaryPart): ...@@ -564,7 +656,7 @@ class Interface(BoundaryPart):
for dof_index in interface_dofs: for dof_index in interface_dofs:
# fem_function.vector() is orderd by dof and not by vertex numbering of the mesh. # 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. # 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()[:]) print("overwritten function \n", fem_function.vector()[:])
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment