diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index d8e99a1bf007ed1e5c2fdf619d64d79f97cf0ca5..0b5678d30b85b631fd5bb2cdbf1e9313b1929fa1 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -162,7 +162,7 @@ class LDDsimulation(object):
'absolute_tolerance': 1E-12,
'relative_tolerance': 1E-10,
'maximum_iterations': 1000,
- 'report': False
+ 'report': True
}
self.debug_solver = debug
@@ -196,6 +196,10 @@ class LDDsimulation(object):
dirichletBC_expression_strings: tp.Dict[int, tp.Dict[str, str]],#
write2file: tp.Dict[str, bool],#
exact_solution: tp.Dict[int, tp.Dict[str, str]] = None, #
+ densities: tp.Dict[int, tp.Dict[str, float]] = None,#
+ include_gravity: bool = False,#
+ gravity_acceleration: float = 9.81,
+
)-> None:
""" set parameters of an instance of class LDDsimulation"""
@@ -224,6 +228,12 @@ class LDDsimulation(object):
self.dirichletBC_expression_strings = dirichletBC_expression_strings
self.exact_solution = exact_solution
self.write2file = write2file
+ # flag to toggle wether or not to include the gravity term in the calculation.
+ self.include_gravity = include_gravity
+ self.gravity_acceleration = gravity_acceleration
+ self.densities = densities
+ if self.include_gravity and self.densities is None:
+ raise RuntimeError("include_gravity = True but no densities were given.")
self._parameters_set = True
def initialise(self) -> None:
@@ -640,7 +650,8 @@ class LDDsimulation(object):
Lam = subdomain.lambda_param['wetting']
name1 ="subdomain" + "{}".format(subdom_ind)
# only show three digits of the mesh size.
- name2 = "_dt" + "{}".format(tau)+"_hmax" + "{number:.{digits}f}".format(number=h, digits=3)
+ name2 = "_dt" + "{}".format(tau) +"_hmax" + "{number:.{digits}f}".format(number=h, digits=3) \
+ + "_gravitiy_" + "{}".format(self.include_gravity)
name3 = "_lambda" + "{}".format(Lam) + "_Lw" + "{}".format(L["wetting"])
filename_param_part = name1 + name2 + name3
filename = self.output_dir + filename_param_part + "_solution" +".xdmf"
@@ -650,6 +661,7 @@ class LDDsimulation(object):
filename_param_part = "subdomain" + "{}".format(subdom_ind)\
+"_dt" + "{}".format(tau)\
+"_hmax" + "{number:.{digits}f}".format(number=h, digits=3)\
+ + "_gravitiy_" + "{}".format(self.include_gravity)\
+"_lambda_w" + "{}".format(Lam_w)\
+"_lambda_nw" + "{}".format(Lam_nw)\
+"_Lw" + "{}".format(L["wetting"])\
@@ -903,6 +915,10 @@ class LDDsimulation(object):
for subdom_num, isR in self.isRichards.items():
interface_list = self._get_interfaces(subdom_num)
# print(f"has_interface list for subdomain {subdom_num}:", interface_list)
+ if self.densities is None:
+ subdom_densities = None
+ else:
+ subdom_densities = self.densities[subdom_num]
self.subdomain.update(#
{subdom_num : dp.DomainPatch(#
subdomain_index = subdom_num,#
@@ -913,11 +929,15 @@ class LDDsimulation(object):
outer_boundary_def_points = self.outer_boundary_def_points[subdom_num],#
interfaces = self.interface,#
has_interface = interface_list,#
+ densities = subdom_densities,#
+ include_gravity = self.include_gravity,#
+ gravity_acceleration = self.gravity_acceleration,#
L = self.L[subdom_num],#
lambda_param = self.lambda_param[subdom_num],#
relative_permeability = self.relative_permeability[subdom_num],#
saturation = self.saturation[subdom_num],#
timestep_size = self.timestep_size,#
+ interpolation_degree = self.interpolation_degree,#
tol = self.tol)
}
)
@@ -1110,8 +1130,11 @@ class LDDsimulation(object):
""" function to create a new output directory for given set of input
parameters and adjust self.output_dir accordingly
"""
+ gravity_string = ""
+ if self.include_gravity:
+ gravity_string = "_with_gravity"
dirname = "mesh_res" + "{}_".format(self.mesh_resolution)\
- + "dt" + "{}".format(self.timestep_size) + "/"
+ + "dt" + "{}".format(self.timestep_size)+ gravity_string + "/"
self.output_dir += dirname
if not os.path.exists(self.output_dir):
diff --git a/LDDsimulation/domainPatch.py b/LDDsimulation/domainPatch.py
index 603826505dffb7513604a5806fc4247df906a4fb..c0aad7b3c9d5ded1981304e9ae6312e6c2289914 100644
--- a/LDDsimulation/domainPatch.py
+++ b/LDDsimulation/domainPatch.py
@@ -75,6 +75,13 @@ class DomainPatch(df.SubDomain):
two-phase flow is assumed on the
patch) we have two parameters, one
for each equation.
+ densities tp.Dict[int, tp.Dict[str, float]] Dictionary for both phases
+ holding the densities for the gravity
+ term.
+ include_gravity bool flag to toggle wether or not to
+ include the gravity term
+ gravity_acceleration float gravity acceleration for the gravity
+ term
## Public Variables
@@ -121,6 +128,10 @@ class DomainPatch(df.SubDomain):
relative_permeability: tp.Dict[str, tp.Callable[...,None]],#
saturation: tp.Callable[..., None],#
timestep_size: float,#
+ densities: tp.Dict[str, float] = None,
+ include_gravity: bool = False,#
+ gravity_acceleration: float = 9.81,
+ interpolation_degree: int = 10,
tol = None,#
):
# because we declare our own __init__ method for the derived class BoundaryPart,
@@ -142,13 +153,16 @@ class DomainPatch(df.SubDomain):
self.outer_boundary_def_points = outer_boundary_def_points
self.interface = interfaces
self.has_interface = has_interface
+ self.densities = densities
+ self.include_gravity = include_gravity
+ self.gravity_acceleration = gravity_acceleration
self.L = L
self.lambda_param = lambda_param
self.relative_permeability = relative_permeability
self.saturation = saturation
# timestep size, tau in the paper
self.timestep_size = timestep_size
-
+ self.interpolation_degree = interpolation_degree
### Class variables set by methods
# sometimes we need to loop over the phases and the length of the loop is
# determined by the model we are assuming
@@ -197,6 +211,11 @@ class DomainPatch(df.SubDomain):
self._init_measures()
self._calc_dof_indices_of_interfaces()
self._calc_interface_has_common_dof_indices()
+ if self.include_gravity:
+ self.gravity_term = dict()
+ self._calc_gravity_expressions()
+ else:
+ self.gravity_term = None
# Dictionary of Arrays (one for 'wetting' and one for 'nonwetting')
# corresponding to the dofs of an FEM function which is initiated
# by calc_gli_term() holding the dofs of the assembled dt*gli*v*ds terms as
@@ -277,6 +296,8 @@ class DomainPatch(df.SubDomain):
S = self.saturation
mu_a = self.viscosity[phase]
source_a = self.source[phase]
+ if self.include_gravity:
+ z_a = self.gravity_term[phase]
# the first part of the form and rhs are the same for both wetting and nonwetting
form1 = (La*pa*phi_a)*dx
@@ -296,6 +317,10 @@ class DomainPatch(df.SubDomain):
# form2 is different for wetting and nonwetting
form2 = dt/mu_a*df.dot(ka(S(pc_prev_iter))*df.grad(pa), df.grad(phi_a))*dx
rhs2 = -(porosity*(S(pc_prev_iter) - S(pc_prev_timestep))*phi_a)*dx
+ if self.include_gravity:
+ # z_a included the minus sign already, see self._calc_gravity_expressions
+ rhs_gravity = -dt/mu_a*df.dot(ka(S(pc_prev_iter))*df.grad(z_a), df.grad(phi_a))*dx
+
elif phase == "nonwetting":
# gather interface forms
interface_forms = []
@@ -314,12 +339,20 @@ class DomainPatch(df.SubDomain):
form2 = dt/mu_a*df.dot(ka(1 - S(pc_prev_iter))*df.grad(pa), df.grad(phi_a))*dx
# the non wetting phase has a + before instead of a - before the bracket
rhs2 = (porosity*(S(pc_prev_iter) - S(pc_prev_timestep))*phi_a)*dx
+ if self.include_gravity:
+ # z_a included the minus sign already, see self._calc_gravity_expressions
+ rhs_gravity = -dt/mu_a*df.dot(ka(1 - S(pc_prev_iter))*df.grad(z_a), df.grad(phi_a))*dx
else:
raise RuntimeError('missmatch for input parameter phase.',
'Expected either phase == wetting or phase == nonwetting ')
return None
form = form1 + form2 + sum(interface_forms)
- rhs_without_gli = rhs1 + rhs2 + dt*source_a*phi_a*dx
+ if self.include_gravity:
+ # z_a included the minus sign already, see self._calc_gravity_expressions
+ rhs_without_gli = rhs1 + rhs2 + dt*source_a*phi_a*dx + rhs_gravity
+ else:
+ rhs_without_gli = rhs1 + rhs2 + dt*source_a*phi_a*dx
+
form_and_rhs = {#
'form': form,#
'rhs_without_gli' : rhs_without_gli
@@ -389,6 +422,7 @@ class DomainPatch(df.SubDomain):
pa_prev = self.pressure_prev_timestep[phase]
# testfunction
phi_a = self.testfunction[phase]
+ z_a = self.gravity_term[phase]
if phase == 'nonwetting' and interface.neighbour_isRichards[subdomain]:
# if the neighbour of our subdomain (with index self.subdomain_index)
# assumes a Richards model, we don't have gli terms for the
@@ -402,12 +436,19 @@ class DomainPatch(df.SubDomain):
# the wetting phase is always present and we always need
# to calculate a gl0 term.
# TODO: add intrinsic permeabilty!!!!
- # TODO: add gravitiy term!!!!!
- flux = -1/mu_a*ka(S(pc_prev))*df.grad(pa_prev)
+ # TODO: add gravity term!!!!!
+ if self.include_gravity:
+ # z_a included the minus sign already, see self._calc_gravity_expressions
+ flux = -1/mu_a*ka(S(pc_prev))*df.grad(pa_prev + z_a)
+ else:
+ flux = -1/mu_a*ka(S(pc_prev))*df.grad(pa_prev)
else:
# phase == 'nonwetting'
# we need anothre flux in this case
- flux = -1/mu_a*ka(1-S(pc_prev))*df.grad(pa_prev)
+ if self.include_gravity:
+ flux = -1/mu_a*ka(1-S(pc_prev))*df.grad(pa_prev + z_a)
+ else:
+ flux = -1/mu_a*ka(1-S(pc_prev))*df.grad(pa_prev)
gl0 = df.dot(flux, n) - Lambda[phase]*pa_prev
gl0_assembled = df.assemble(dt*gl0*phi_a*ds).get_local()
@@ -798,6 +839,17 @@ class DomainPatch(df.SubDomain):
# print(f"python id of self._interface_has_common_dof_indices[{interf_ind}][{phase}]", id(self._interface_has_common_dof_indices[interf_ind][phase]))
+ def _calc_gravity_expressions(self):
+ """ calculate the gravity term if it is included"""
+ g = df.Constant(self.gravity_acceleration) #
+ for phase in self.has_phases:
+ rho = df.Constant(self.densities[phase])
+ self.gravity_term.update(
+ {phase: df.Expression('-g*rho*x[1]', domain = self.mesh, #
+ degree = self.interpolation_degree,
+ g = g, rho = rho)}
+ )
+
def write_solution_to_xdmf(self, file: tp.Type[SolutionFile], #
time: float = None,#
write_iter_for_fixed_time: bool = False,#