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

LDDsimulation.py

parent 9d3d817d
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
LDDsimulation.py
+ 72
− 5
View file @ 54a7c58b
......@@ -84,8 +84,8 @@ class LDDsimulation(object):
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 = []
# dictionary of objects of class DomainPatch initialised by self._init_subdomains()
self.subdomain = dict()
def set_parameters(self,#
output_dir: str,#
......@@ -102,6 +102,8 @@ class LDDsimulation(object):
relative_permeability: tp.Dict[int, tp.Dict[str, tp.Callable[...,None]] ],#
saturation: tp.Dict[int, tp.Callable[...,None]],#
timestep: tp.Dict[int, float],#
sources: tp.Dict[int, tp.Dict[str, str]],#
initial_conditions: tp.Dict[int, tp.Dict[str, str]]
)-> None:
""" set parameters of an instance of class LDDsimulation"""
......@@ -118,6 +120,8 @@ class LDDsimulation(object):
self.relative_permeability = relative_permeability
self.saturation = saturation
self.timestep = timestep
self.sources = sources
self.initial_conditions = initial_conditions
self._parameters_set = True
def initialise(self) -> None:
......@@ -130,7 +134,7 @@ class LDDsimulation(object):
self._init_meshes_and_markers()
self._init_interfaces()
self._init_subdomains()
self._init_initial_values()
## Private methods
def _init_meshes_and_markers(self, subdomain_def_points: tp.List[tp.List[df.Point]] = None,#
......@@ -260,7 +264,8 @@ class LDDsimulation(object):
# 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(#
self.subdomain.update(#
{subdom_num : dp.DomainPatch(#
subdomain_index = subdom_num,#
isRichards = isR,#
mesh = self.mesh_subdomain[subdom_num],#
......@@ -273,7 +278,7 @@ class LDDsimulation(object):
relative_permeability = self.relative_permeability[subdom_num],#
saturation = self.saturation[subdom_num],#
timestep = self.timestep[subdom_num],#
))
)})
......@@ -300,3 +305,65 @@ class LDDsimulation(object):
interface_of_subdomain.append(interface_ind)
return interface_of_subdomain
def _init_initial_values(self, interpolation_degree: int = 2):
""" set initial values
"""
for num, subdomain in self.subdomain.items():
mesh = subdomain.mesh
V = subdomain.function_space
# p0 contains both pw_0 and pnw_0 for subdomain[num]
p0 = self.initial_conditions[num]
if subdomain.isRichards:
# note that the default interpolation degree is 2
pw0 = df.Expression(p0['wetting'], domain = mesh, degree = interpolation_degree)
subdomain.pressure_prev_iter = {'wetting' : df.interpolate(pw0, V['wetting'])}
subdomain.pressure_prev_timestep = {'wetting' : df.interpolate(pw0, V['wetting'])}
else:
pw0 = df.Expression(p0['wetting'], domain = mesh, degree = interpolation_degree)
pnw0 = df.Expression(p0['nonwetting'], domain = mesh, degree = interpolation_degree)
subdomain.pressure_prev_iter = {'wetting' : df.interpolate(pw0, V['wetting']),#
'nonwetting' : df.interpolate(pnw0, V['nonwetting'])}
subdomain.pressure_prev_timestep = {'wetting' : df.interpolate(pw0, V['wetting']),#
'nonwetting' : df.interpolate(pnw0, V['nonwetting'])}
def _eval_sources(self, interpolation_degree: int = 2, time: float = 0):
""" evaluate time dependent source terms or initialise them if time == 0
"""
if np.abs(time) < self.tol:
# here t = 0 and we have to initialise the sources.
for num, subdomain in self.subdomain.items():
mesh = subdomain.mesh
V = subdomain.function_space
# p0 contains both pw_0 and pnw_0 for subdomain[num]
f = self.sources[num]
if subdomain.isRichards:
# note that the default interpolation degree is 2
fw = df.Expression(f['wetting'], domain = mesh, degree = interpolation_degree, t = time)
# subdomain.source = {'wetting' : df.interpolate(fw, V['wetting'])}
subdomain.source = {'wetting' : fw}
else:
fw = df.Expression(f['wetting'], domain = mesh, degree = interpolation_degree, t = time)
fnw = df.Expression(f['nonwetting'], domain = mesh, degree = interpolation_degree, t = time)
# subdomain.pressure_prev_iter = {'wetting' : df.interpolate(fw, V['wetting']),#
# 'nonwetting' : df.interpolate(fnw, V['nonwetting'])}
# subdomain.pressure_prev_timestep = {'wetting' : df.interpolate(fw, V['wetting']),#
# 'nonwetting' : df.interpolate(fnw, V['nonwetting'])}
subdomain.source = {'wetting' : fw,#
'nonwetting' : fnw}
else:
for num, subdomain in self.subdomain.items():
# mesh = subdomain.mesh
# V = subdomain.function_space
# # p0 contains both pw_0 and pnw_0 for subdomain[num]
# f = self.sources[num]
if subdomain.isRichards:
# note that the default interpolation degree is 2
# fw = df.Expression(f['wetting'], domain = mesh, degree = interpolation_degree, t = time)
# subdomain.source = {'wetting' : df.interpolate(fw, V['wetting'])}
subdomain.source['wetting'].t = time
else:
# fw = df.Expression(f['wetting'], domain = mesh, degree = interpolation_degree, t = time)
# fnw = df.Expression(f['nonwetting'], domain = mesh, degree = interpolation_degree, t = time)
subdomain.source['wetting'].t = time
subdomain.source['nonwetting'].t = time
RR-2-patch-test.py
+ 29
− 10
View file @ 54a7c58b
......@@ -115,7 +115,15 @@ sat_pressure_relationship = {#
2: ft.partial(saturation, subdomain_index = 2)
}
source_expression = {
1: {'wetting': '4.0/pow(1 + x[0]*x[0] + x[1]*x[1], 2) - t/sqrt( pow(1 + t*t, 3)*(1 + x[0]*x[0] + x[1]*x[1]) )'},
2: {'wetting': '2.0*(1-x[1]*x[1])/pow(1 + x[1]*x[1], 2) - 2*t/(3*pow( pow(1 + t*t, 4)*(1 + x[1]*x[1]), 1/3))'}
}
initial_condition = {
1: {'wetting': '-(x[0]*x[0] + x[1]*x[1])'},#
2: {'wetting': '-x[1]*x[1]'}
}
# def saturation(pressure, subdomain_index):
# # inverse capillary pressure-saturation-relationship
......@@ -138,7 +146,9 @@ simulation.set_parameters(output_dir = "",#
lambda_param = lambda_param,#
relative_permeability = relative_permeability,#
saturation = sat_pressure_relationship,#
timestep = timestep#
timestep = timestep,#
sources = source_expression,#
initial_conditions = initial_condition,#
)
simulation.initialise()
......@@ -161,6 +171,13 @@ mesh_subdomain = simulation.mesh_subdomain
interface = simulation.interface
interface_marker = simulation.interface_marker
boundary_marker1 = simulation.subdomain[1].boundary_marker
boundary_marker2 = simulation.subdomain[2].boundary_marker
dx_1 = simulation.subdomain[1].dx
dx_2 = simulation.subdomain[2].dx
ds_1 = simulation.subdomain[1].ds
ds_2 = simulation.subdomain[2].ds
subdomain_boundary_marker1 = df.MeshFunction('size_t', mesh_subdomain[1], mesh_subdomain[1].topology().dim()-1)
subdomain_boundary_marker2 = df.MeshFunction('size_t', mesh_subdomain[2], mesh_subdomain[2].topology().dim()-1)
subdomain_boundary_marker1.set_all(0)
......@@ -169,9 +186,11 @@ subdomain_boundary_marker2.set_all(0)
interface[0].mark(subdomain_boundary_marker1, 1)
interface[0].mark(subdomain_boundary_marker2, 1)
# subdomain_boundary_marker1 = simulation.subdomain[1].boundary_marker
# domain measures
# dx1 = simulation.subdomain[1].dx
dx1 = df.Measure('dx', domain = mesh_subdomain[1], subdomain_data = subdomain_boundary_marker1)
dx2 = df.Measure('dx', domain = mesh_subdomain[2], subdomain_data = subdomain_boundary_marker2)
......@@ -275,17 +294,15 @@ interface[0].write_dofs(from_function = f1test, #
interface_dofs = dofs_on_interface1,#
dof_to_vert_map = dom1_d2v,#
local_to_parent_vertex_map = dom1_parent_vertex_indices,#
phase = 'wetting',#
pressure = True,#
phase = 'wetting',
subdomain_ind = 1)
interface[0].read_dofs(to_function = f2test, #
interface_dofs = dofs_on_interface2,#
dof_to_vert_map = dom2_d2v,#
local_to_parent_vertex_map = dom2_parent_vertex_indices,#
phase = 'wetting',#
pressure = True,#
subdomain_ind = 1)
phase = 'wetting',
subdomain_ind = 2)
i=0
for x in coordinates2:
......@@ -342,8 +359,11 @@ outerBC2 = df.DirichletBC(V2, p2_exact, subdomain_boundary_marker2, 0)
### source terms
source1 = df.Expression('4.0/pow(1 + x[0]*x[0] + x[1]*x[1], 2) - t/sqrt( pow(1 + t*t, 3)*(1 + x[0]*x[0] + x[1]*x[1]) )', domain = mesh_subdomain[1], degree = 2, t = 0)
# source1 = df.Expression('4.0/pow(1 + x[0]*x[0] + x[1]*x[1], 2) - t/sqrt( pow(1 + t*t, 3)*(1 + x[0]*x[0] + x[1]*x[1]) )', domain = mesh_subdomain[1], degree = 2, t = 0)
# source2 = df.Expression('2.0*(1-x[1]*x[1])/pow(1 + x[1]*x[1], 2) - 2*t/(3*pow( pow(1 + t*t, 4)*(1 + x[1]*x[1]), 1/3))', domain = mesh_subdomain[2], degree = 2, t = 0)
source1 = df.Expression(source_expression[1]['wetting'], domain = mesh_subdomain[1], degree = 2, t = 0)
source2 = df.Expression('2.0*(1-x[1]*x[1])/pow(1 + x[1]*x[1], 2) - 2*t/(3*pow( pow(1 + t*t, 4)*(1 + x[1]*x[1]), 1/3))', domain = mesh_subdomain[2], degree = 2, t = 0)
source1h = df.interpolate(source1, V1)
source2h = df.interpolate(source2, V2)
......@@ -361,7 +381,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*v1)*dx1 -(timestep*(g1_i)*v1)*ds1
rhs1 = (L1*p1_i*v1)*dx1 - ((saturation(p1_i, 1) - saturation(p1_0, 1))*v1)*dx1 + (timestep*(source1 - g1_i)*v1)*ds1
A1 = df.assemble(a1)
b1 = df.assemble(rhs1)
p_gamma1.apply(A1,b1)
......@@ -376,7 +396,6 @@ for iterations in range(1):
dof_to_vert_map = dom1_d2v,#
local_to_parent_vertex_map = dom1_parent_vertex_indices,#
phase = 'wetting',#
pressure = True,#
subdomain_ind = 1)
#
# Save mesh to file
......
domainPatch.py
+ 11
− 9
View file @ 54a7c58b
......@@ -143,10 +143,11 @@ class DomainPatch(df.SubDomain):
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
self.source = None
# dictionary holding the previous iteration.
self.pressure_prev_iter = None
# dictionary holding the pressures of the previous timestep t_n.
self.pressure_prev_timestep = None
# test function(s) for the form set by _init_function_space
self.TestFunction = None
self._init_function_space()
......@@ -168,8 +169,8 @@ class DomainPatch(df.SubDomain):
# 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']
pw_prev_iter = self.pressure_prev_iter['wetting']
pw_prev_timestep = self.pressure_prev_timestep['wetting']
v = self.testfunction['wetting']
Lambda = self.lambda_param['wetting']
kw = self.relative_permeability['wetting']
......@@ -552,8 +553,9 @@ class Interface(BoundaryPart):
def init_interface_values(self,#
interface_marker: df.MeshFunction,#
interface_marker_value: int) -> dict:
""" allocate dictionary that is used to store the interface dof values
as pairs(parent_mesh_index: dof_value)
""" allocate dictionaries that are used to store the interface dof values
of the pressure, previous pressure, gl decoupling term as well as previous
gl update term respectively as pairs(parent_mesh_index: dof_value)
The MeshFunction interface_marker is supposed to be one defined on the
global mesh, so that the indices used in the dictionary are the global
......@@ -595,6 +597,7 @@ class Interface(BoundaryPart):
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})
# end init_interface_values
def write_dofs(self, from_function: df.Function, #
interface_dofs: np.array,#
......@@ -784,7 +787,6 @@ class Interface(BoundaryPart):
if self._is_on_boundary_part(x):
# print(f"Vertex {x} with index {vert_num} is on interface")
# print(f"interface_vertex_number = {interface_vertex_number}")
vertex_indices[interface_vertex_number] = vert_num
interface_vertex_number += 1
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment