diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index c6f9da6d1afe00104434ee783f868567fe9d6270..0fa2541099e744aaba81deb800d2476c11656b81 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -81,9 +81,9 @@ class LDDsimulation(object):
## Private variables
# maximal number of L-iterations that the LDD solver uses.
- self._max_iter_num = 100
+ self._max_iter_num = 5
# TODO rewrite this with regard to the mesh sizes
- self.calc_tol = 10^-6
+ self.calc_tol = self.tol
# 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
@@ -111,6 +111,29 @@ class LDDsimulation(object):
# These are interpolated and stored to each respective subdomain by
# self._init_initial_values()
self.initial_conditions = None
+ ########## SOLVER DEFINITION AND PARAMETERS ########
+ # print("\nLinear Algebra Backends:")
+ # df.list_linear_algebra_backends()
+ # print("\nLinear Solver Methods:")
+ # df.list_linear_solver_methods()
+ # print("\nPeconditioners for Krylov Solvers:")
+ # df.list_krylov_solver_preconditioners()
+
+ ### Define the linear solver to be used.
+ solver = 'bicgstab' # biconjugate gradient stabilized method
+ preconditioner = 'ilu' # incomplete LU factorization
+ self.linear_solver = df.KrylovSolver(solver, preconditioner)
+ # we use the previous iteration as an initial guess for the linear solver.
+ solver_param = self.linear_solver.parameters
+ solver_param['nonzero_initial_guess'] = True
+ # solver_param['absolute_tolerance'] = 1E-12
+ # solver_param['relative_tolerance'] = 1E-9
+ solver_param['maximum_iterations'] = 1000
+ solver_param['report'] = True
+ # ## print out set solver parameters
+ # for parameter, value in self.linear_solver.parameters.items():
+ # print(f"parameter: {parameter} = {value}")
+ df.info(solver_param, True)
def set_parameters(self,#
output_dir: str,#
@@ -181,7 +204,7 @@ class LDDsimulation(object):
# before the iteration gets started all iteration numbers must be nulled
# and the subdomain_calc_isdone dictionary must be set to False
for ind, subdomain in self.subdomain.items():
- subdomain_calc_isdone.update({ind, False})
+ subdomain_calc_isdone.update({ind: False})
# reset all interface[has_interface].current_iteration[ind] = 0, for
# index has_interface in subdomain.has_interface
subdomain.null_all_interface_iteration_numbers()
@@ -194,15 +217,15 @@ class LDDsimulation(object):
# the following only needs to be done when time is not equal 0 since
# our initialisation functions already took care of setting what follows
# for the initial iteration step at time = 0.
- if time > self.tol:
+ if time > self.calc_tol:
# this is the beginning of the solver, so iteration must be set
# to 0.
subdomain.iteration_number = 0
- # set the solution of the previous timestep as new prev_timestep pressure
- subdomain.pressure_prev_timestep = subdomain.pressure
- # TODO recheck if
- subdomain.pressure_prev_iter = subdomain.pressure
- # needs to be set also
+ for phase in subdomain.has_phases:
+ # set the solution of the previous timestep as new prev_timestep pressure
+ subdomain.pressure_prev_timestep[phase].vector().set_local(
+ subdomain.pressure[phase].vector().get_local()
+ )
### actual iteration starts here
# gobal stopping criterion for the iteration.
@@ -220,21 +243,24 @@ class LDDsimulation(object):
subdomain.iteration_number = iteration
# solve the problem on subdomain
self.Lsolver_step(subdomain_index = sd_index,#
- debug = True
- )
- subsequent_iterations_err = self.calc_iteration_error(
- subdomain_index = sd_index,#
- error_type = "L2"
+ debug = True,
)
- # check if error criterion has been met
- if subsequent_iterations_err < self.calc_tol:
- subdomain_calc_isdone[sd_index] = True
+ # subsequent_iterations_err = self.calc_iteration_error(
+ # subdomain_index = sd_index,#
+ # error_type = "L2"
+ # )
+ # # check if error criterion has been met
+ # if subsequent_iterations_err < self.calc_tol:
+ # subdomain_calc_isdone[sd_index] = True
# prepare next iteration
# write the newly calculated solution to the inteface dictionaries
# for communication
subdomain.write_pressure_to_interfaces()
- subdomain.pressure_prev_iter = subdomain.pressure
+ for phase in subdomain.has_phases:
+ subdomain.pressure_prev_iter[phase].vector().set_local(
+ subdomain.pressure[phase].vector().get_local()
+ )
else:
# one subdomain is done. Check if all subdomains are done to
# stop the while loop if needed.
@@ -255,7 +281,7 @@ class LDDsimulation(object):
iteration += 1
# end iteration while loop.
- def Lsolver_step(self, subdomain_index: int],#
+ def Lsolver_step(self, subdomain_index: int,#
debug: bool = False) -> None:
""" L-scheme solver iteration step for an object of class subdomain
@@ -273,28 +299,54 @@ class LDDsimulation(object):
"""
subdomain = self.subdomain[subdomain_index]
iteration = subdomain.iteration_number
- if subdomain.isRichards:
+ # calculate the gli terms on the right hand side and store
+ subdomain.calc_gli_term()
+
+ for phase in subdomain.has_phases:
# extract L-scheme form and rhs (without gli term) from subdomain.
- governing_problem = subdomain.governing_problem(phase = 'wetting')
- form = governing_problem[form]
- rhs_without_gli = governing_problem[rhs_without_gli]
+ governing_problem = subdomain.governing_problem(phase = phase)
+ form = governing_problem['form']
+ rhs_without_gli = governing_problem['rhs_without_gli']
# assemble the form and rhs
form_assembled = df.assemble(form)
rhs_without_gli_assembled = df.assemble(rhs_without_gli)
- # calculate the gli term
- gli_term_assembled = subdomain.calc_gli_term(iteration = iteration)
+ # access the assembled gli term for the rhs.
+ gli_term_assembled = subdomain.gli_assembled[phase]
# subdomain.calc_gli_term() asslembles gli but on the left hand side
# so gli_term_assembled needs to actually be subtracted from the rhs.
- rhs_assembled = rhs_without_gli_assembled - gli_term_assembled
+ if debug and subdomain.mesh.num_cells() < 36:
+ print("\nSystem before applying outer boundary conditions:")
+ print(f"phase = {phase}: rhs_without_gli:\n", rhs_without_gli_assembled.get_local())
+ rhs_without_gli_assembled.set_local(rhs_without_gli_assembled.get_local() - gli_term_assembled)
+ rhs_assembled = rhs_without_gli_assembled
+
+ if debug and subdomain.mesh.num_cells() < 36:
+ # print(f"phase = {phase}: form_assembled:\n", form_assembled.array())
+ # print(f"phase = {phase}: rhs_without_gli:\n", rhs_without_gli_assembled.get_local())
+ print(f"phase = {phase}: gli_term:\n", gli_term_assembled)
+ print(f"phase = {phase}: rhs_assembled:\n", rhs_assembled.get_local())
+
# apply outer Dirichlet boundary conditions if present.
if subdomain.outer_boundary is not None:
- self.outerBC[subdomain_index]['wetting'].apply(form_assembled, rhs_assembled)
- solution_pw = subdomain.pressure['wetting'].vector()
+ self.outerBC[subdomain_index][phase].apply(form_assembled, rhs_assembled)
- else:
- print("implement two phase assembly")
+ if debug and subdomain.mesh.num_cells() < 36:
+ print("\nSystem after applying outer boundary conditions:")
+ # print(f"phase = {phase}: form_assembled:\n", form_assembled.array())
+ print(f"phase = {phase}: rhs_assembled:\n", rhs_assembled.get_local())
+ # # set previous iteration as initial guess for the linear solver
+ # pi_prev = subdomain.pressure_prev_iter[phase].vector().get_local()
+ # subdomain.pressure[phase].vector().set_local(pi_prev)
+ if debug and subdomain.mesh.num_cells() < 36:
+ print("\npressure before solver:\n", subdomain.pressure[phase].vector().get_local())
+
+ solver = self.linear_solver
+ solver.solve(form_assembled, subdomain.pressure[phase].vector(), rhs_assembled)
+
+ if debug and subdomain.mesh.num_cells() < 36:
+ print("\npressure after solver:\n", subdomain.pressure[phase].vector().get_local())
## Private methods
def _init_meshes_and_markers(self, subdomain_def_points: tp.List[tp.List[df.Point]] = None,#
@@ -358,7 +410,7 @@ class LDDsimulation(object):
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
+ """ generate interfaces and interface markerfunctions for suddomains and
set the internal lists used by the LDDsimulation class.
This initialises a list of interfaces and global interface marker functions
@@ -405,7 +457,7 @@ class LDDsimulation(object):
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,#
+ tol=self.tol,#mesh.hmin()/100,
internal=True,#
adjacent_subdomains = adjacent_subdomains[num],#
isRichards = self.isRichards))#
@@ -423,6 +475,7 @@ class LDDsimulation(object):
# Therefor it is used in the subdomain initialisation loop.
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)
self.subdomain.update(#
{subdom_num : dp.DomainPatch(#
subdomain_index = subdom_num,#
@@ -438,6 +491,7 @@ class LDDsimulation(object):
relative_permeability = self.relative_permeability[subdom_num],#
saturation = self.saturation[subdom_num],#
timestep_size = self.timestep_size[subdom_num],#
+ tol = self.tol
)})
@@ -450,7 +504,7 @@ class LDDsimulation(object):
This method determins all (global) indices of interfaces belonging to
subdomain with index subdomain_index from adjacent_subdomains.
"""
- if not adjacent_subdomains:
+ if adjacent_subdomains is None:
adjacent_subdomains = self.adjacent_subdomains
else:
# overwrite what has been set by init()
@@ -480,24 +534,45 @@ class LDDsimulation(object):
if subdomain.isRichards:
# note that the default interpolation degree is 2
pw0 = df.Expression(p0['wetting'], domain = mesh, degree = interpolation_degree)
+ subdomain.pressure = {
+ 'wetting': df.Function(V['wetting'])
+ }
+ subdomain.pressure_prev_timestep = {
+ 'wetting': df.Function(V['wetting'])
+ }
subdomain.pressure_prev_iter = {'wetting' : df.interpolate(pw0, V['wetting'])}
- subdomain.pressure_prev_timestep = {'wetting' : df.interpolate(pw0, V['wetting'])}
+ # print("vector()",subdomain.pressure_prev_iter['wetting'].vector().get_local())
+ pw_prev = subdomain.pressure_prev_iter['wetting'].vector().get_local()
+ subdomain.pressure_prev_timestep['wetting'].vector().set_local(pw_prev)
+ subdomain.pressure['wetting'].vector().set_local(pw_prev)
else:
pw0 = df.Expression(p0['wetting'], domain = mesh, degree = interpolation_degree)
pnw0 = df.Expression(p0['nonwetting'], domain = mesh, degree = interpolation_degree)
+ subdomain.pressure = {
+ 'wetting': df.Function(V['wetting']),#
+ 'nonwetting': df.Function(V['nonwetting'])
+ }
+ subdomain.pressure_prev_timestep = {
+ 'wetting': df.Function(V['wetting']),
+ 'nonwetting': df.Function(V['nonwetting'])
+ }
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'])}
+
+ pw_prev = subdomain.pressure_prev_iter['wetting'].vector().get_local()
+ pnw_prev = subdomain.pressure_prev_iter['nonwetting'].vector().get_local()
+ subdomain.pressure_prev_timestep['wetting'].vector().set_local(pw_prev)
+ subdomain.pressure['wetting'].vector().set_local(pw_prev)
+ subdomain.pressure_prev_timestep['nonwetting'].vector().set_local(pnw_prev)
+ subdomain.pressure['nonwetting'].vector().set_local(pnw_prev)
# populate interface dictionaries with pressure values.
subdomain.write_pressure_to_interfaces(initial_values = True)
- def update_DirichletBC_dictionary(self,
+ def update_DirichletBC_dictionary(self, subdomain_index: int,#
interpolation_degree: int = 2, #
time: float = 0,#
- # exact_solution: bool = False,#
- subdomain_index: int):
+ ):
""" update time of the dirichlet boundary object for subdomain with
index subdomain_index.
@@ -512,7 +587,7 @@ class LDDsimulation(object):
mesh = subdomain.mesh
V = subdomain.function_space
# Here the dictionary has to be created first because t = 0.
- if np.abs(time) < self.tol :
+ if np.abs(time) < self.calc_tol:
self.outerBC.update({num: dict()})
self.dirichletBC_dfExpression.update({num: dict()})
@@ -529,7 +604,7 @@ class LDDsimulation(object):
self.outerBC[num].update({'wetting': df.DirichletBC(V['wetting'], pDCw, boundary_marker, 1)})
else:
- if np.abs(time) < self.tol :
+ if np.abs(time) < self.calc_tol :
# time = 0 so the Dolfin Expression has to be created first.
pDCw = df.Expression(pDC['wetting'], domain = mesh, degree = interpolation_degree, t=time)
pDCnw = df.Expression(pDC['nonwetting'], domain = mesh, degree = interpolation_degree, t=time)
@@ -550,13 +625,14 @@ class LDDsimulation(object):
print("subdomain is an inner subdomain and has no outer boundary.\n",
"no Dirichlet boundary has been set.")
- def _eval_sources(self, interpolation_degree: int = 2, #
- time: float = 0,
- subdomain_index: int):
+ def _eval_sources(self, subdomain_index: int, #
+ interpolation_degree: int = 2, #
+ time: float = 0,#
+ ):
""" evaluate time dependent source terms or initialise them if time == 0
"""
subdomain = self.subdomain[subdomain_index]
- if np.abs(time) < self.tol:
+ if np.abs(time) < self.calc_tol:
# here t = 0 and we have to initialise the sources.
mesh = subdomain.mesh
V = subdomain.function_space
diff --git a/LDDsimulation/domainPatch.py b/LDDsimulation/domainPatch.py
index 1aad5c959b2df321a2b9bec2a65e8009d8149de0..f318f15689514e52ff0b3d21a4facbbc0a4aa587 100644
--- a/LDDsimulation/domainPatch.py
+++ b/LDDsimulation/domainPatch.py
@@ -117,12 +117,18 @@ class DomainPatch(df.SubDomain):
relative_permeability: tp.Dict[str, tp.Callable[...,None]],#
saturation: tp.Callable[..., None],#
timestep_size: float,#
+ tol = None,#
):
# because we declare our own __init__ method for the derived class BoundaryPart,
# we overwrite the __init__-method of the parent class df.SubDomain. However,
# we need the methods from the parent class, so we call the __init__-method
# from df.SubDomain to access it.
df.SubDomain.__init__(self)
+ if tol:
+ self.tol = tol
+ else:
+ self.tol = df.DOLFIN_EPS
+
### Class Variables/Objects set by the input to the constructor
self.subdomain_index = subdomain_index
self.isRichards = isRichards
@@ -168,12 +174,19 @@ class DomainPatch(df.SubDomain):
# a sparse vector (here, dt = self.timestep_size). This vector is read by the solver and added to the rhs
# which is given by self.govering_problem().
#
- self.gli_assembled = None
+ self.gli_assembled = dict()
# Dictionary of FEM function of self.function_space holding the interface
# dofs of the previous iteration of the neighbouring pressure. This is used
# to assemble the gli terms in the method self.calc_gli_term().
self.neighbouring_interface_pressure = None
+ # sometimes we need to loop over the phases and the length of the loop is
+ # determined by the model we are assuming
+ if self.isRichards:
+ self.has_phases =['wetting']
+ else:
+ self.has_phases =['wetting', 'nonwetting']
+
self._init_function_space()
self._init_dof_and_vertex_maps()
self._init_markers()
@@ -202,13 +215,8 @@ class DomainPatch(df.SubDomain):
else:
from_func = self.pressure
- if self.isRichards:
- subdomain_has_phases = ['wetting']
- else:
- subdomain_has_phases = ['wetting', 'nonwetting']
-
for ind in self.has_interface:
- for phase in subdomain_has_phases:
+ for phase in self.has_phases:
self.interface[ind].read_pressure_dofs(from_function = from_func[phase], #
interface_dofs = self._dof_indices_of_interface[ind][phase],#
dof_to_vert_map = self.dof2vertex[phase],#
@@ -217,7 +225,7 @@ class DomainPatch(df.SubDomain):
subdomain_ind = self.subdomain_index)
- def govering_problem(self, phase: str) -> tp.Dict[str, fl.Form]:
+ def governing_problem(self, phase: str) -> tp.Dict[str, fl.Form]:
""" return the governing form and right hand side for phase phase as a dictionary.
return the governing form ant right hand side for phase phase (without the gli terms) as a dictionary.
@@ -230,108 +238,76 @@ class DomainPatch(df.SubDomain):
dt = self.timestep_size
porosity = self.porosity
if self.isRichards:
+ # this assumes that atmospheric pressure has been normalised to 0.
+ pc_prev_iter = self.pressure_prev_iter['wetting']
+ pc_prev_timestep = self.pressure_prev_timestep['wetting']
if phase == 'nonwetting':
print("CAUTION: You invoked domainPatch.governing_problem() with\n", #
"Parameter phase = 'nonwetting'. But isRichards = True for",
- "current subdomain. \nReturning wetting phase equation only.\n")
- 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.pressure_prev_iter['wetting']
- pw_prev_timestep = self.pressure_prev_timestep['wetting']
- phi_w = self.testfunction['wetting']
- Lambda = self.lambda_param['wetting']
- kw = self.relative_permeability['wetting']
- S = self.saturation
- mu_w = self.viscosity['wetting']
- source_w = self.source['wetting']
- # we need to have all interfaces in the form
- interface_forms = []
- for interface in self.has_interface:
- interface_forms.append((dt*Lambda*pw*phi_w)*ds(interface))
- form1 = Lw*pw*phi_w*dx
- form2 = dt/mu_w*df.dot(kw(S(pw_prev_iter))*df.grad(pw), df.grad(v))*dx
- form = form1 + form2 + df.sum(interface_forms)
- # # assemble rhs
- rhs1 = (Lw*pw_prev_iter*phi_w)*dx
- rhs2 = -(porosity*(S(pw_prev_iter) - S(pw_prev_timestep))*phi_w)*dx
- rhs_without_gli = rhs1 + rhs2 + dt*source_w*phi_w*dx
- form_and_rhs = {#
- 'form': form,#
- 'rhs_without_gli' : rhs_without_gli
- }
- return form_and_rhs
+ "current subdomain. \nResetting phase = 'wetting'.\n")
+ phase = 'wetting'
else:
- La = self.L['phase']
- pa = self.pressure['phase']
- # this is pw_{i-1} in the paper
pw_prev_iter = self.pressure_prev_iter['wetting']
pnw_prev_iter = self.pressure_prev_iter['nonwetting']
pw_prev_timestep = self.pressure_prev_timestep['wetting']
pnw_prev_timestep = self.pressure_prev_timestep['nonwetting']
pc_prev_iter = pnw_prev_iter - pw_prev_iter
pc_prev_timestep = pnw_prev_timestep - pw_prev_timestep
- phi_a = self.testfunction['phase']
- Lambda_a = self.lambda_param['phase']
- ka = self.relative_permeability['phase']
- S = self.saturation
- mu_a = self.viscosity['phase']
- source_a = self.source['phase']
- # the forms of wetting and nonwetting phase differ by a minus sign
- # and by interface_form terms if the neighbour of the current patch
- # assumes Richards model
- if phase == "wetting":
- # gather interface forms
- interface_forms = []
- for interface in self.has_interface:
- interface_forms.append((dt*Lambda_a*pa*phi_a)*ds(interface))
-
- form1 = (La*pa*phi_a)*dx
- form2 = dt/mu_a*df.dot(ka(S(pc_prev_iter))*df.grad(pa), df.grad(phi_a))*dx
- form = form1 + form2 + df.sum(interface_forms)
-
- rhs1 = (La*pw_prev_iter*phi_a)*dx
- rhs2 = -(porosity*(S(pc_prev_iter) - S(pc_prev_timestep))*phi_a)*dx
- rhs_without_gli = rhs1 + rhs2 + dt*source_a*phi_a*dx
-
- form_and_rhs = {#
- 'form': form,#
- 'rhs_without_gli' : rhs_without_gli
- }
- return form_and_rhs
-
- elif phase == "nonwetting":
- # gather interface forms
- interface_forms = []
- for interface in self.has_interface:
- if self.interface[interface].neighbour_isRichards[self.subdomain_index]:
- # if the neighbour of our subdomain (with index self.subdomain_index)
- # assumes a Richards model, we assume constant normalised atmospheric pressure
- # and no interface term is practically appearing.
- interface_forms.append((df.Constant(0)*phi_a)*ds(interface))
- else:
- interface_forms.append((dt*Lambda_a*pa*phi_a)*ds(interface))
-
- form1 = (La*pa*phi_a)*dx
- form2 = dt/mu_a*df.dot(ka(1 - S(pc_prev_iter))*df.grad(pa), df.grad(phi_a))*dx
- form = form1 + form2 + df.sum(interface_forms)
-
- rhs1 = (La*pnw_prev_iter*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
- rhs_without_gli = rhs1 + rhs2 + dt*source_a*phi_a*dx
- form_and_rhs = {#
- 'form': form,#
- 'rhs_without_gli' : rhs_without_gli
- }
- return form_and_rhs
+ La = self.L[phase]
+ pa = self.trialfunction[phase]
+ # this is pw_{i-1} in the paper
+ pa_prev_iter = self.pressure_prev_iter[phase]
+ phi_a = self.testfunction[phase]
+ Lambda_a = self.lambda_param[phase]
+ ka = self.relative_permeability[phase]
+ S = self.saturation
+ mu_a = self.viscosity[phase]
+ source_a = self.source[phase]
+
+ # the first part of the form and rhs are the same for both wetting and nonwetting
+ form1 = (La*pa*phi_a)*dx
+ rhs1 = (La*pa_prev_iter*phi_a)*dx
+ # the forms of wetting and nonwetting phase differ by a minus sign
+ # and by interface_form terms if the neighbour of the current patch
+ # assumes Richards model
+ if phase == "wetting":
+ # gather interface forms
+ interface_forms = []
+ for interface in self.has_interface:
+ # interface_forms += (dt*Lambda_a*pa*phi_a)*ds(interface)
+ interface_forms.append((dt*Lambda_a*pa*phi_a)*ds(interface))
+ # 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
+ elif phase == "nonwetting":
+ # gather interface forms
+ interface_forms = []
+ for interface in self.has_interface:
+ if self.interface[interface].neighbour_isRichards[self.subdomain_index]:
+ # if the neighbour of our subdomain (with index self.subdomain_index)
+ # assumes a Richards model, we assume constant normalised atmospheric pressure
+ # and no interface term is practically appearing.
+ # interface_forms += (df.Constant(0)*phi_a)*ds(interface)
+ interface_forms.append((df.Constant(0)*phi_a)*ds(interface))
+ else:
+ # interface_forms += (dt*Lambda_a*pa*phi_a)*ds(interface)
+ interface_forms.append((dt*Lambda_a*pa*phi_a)*ds(interface))
+ 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
+ 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
+ form_and_rhs = {#
+ 'form': form,#
+ 'rhs_without_gli' : rhs_without_gli
+ }
+ return form_and_rhs
- else:
- raise RuntimeError('missmatch for input parameter phase.',
- 'Expected either phase == wetting or phase == nonwetting ')
- return None
def calc_gli_term(self) -> None: # tp.Dict[str, fl.Form]
"""calculate the gl terms for the iteration number of the subdomain
@@ -354,13 +330,14 @@ class DomainPatch(df.SubDomain):
# into one array.
gli_assembled_tmp = dict()
if self.isRichards:
+ # print("number of dofs: ", self.pressure_prev_iter['wetting'].vector()[:].size)
gli_assembled_tmp = {
- 'wetting': np.zeros(subdomain.pressure['wetting'].vector().size, dtype=float)
+ 'wetting': np.zeros(self.pressure['wetting'].vector()[:].size, dtype=float)
}
else:
gli_assembled_tmp = {
- 'wetting': np.zeros(subdomain.pressure['wetting'].vector().size, dtype=float),
- 'nonwetting': np.zeros(subdomain.pressure['nonwetting'].vector().size, dtype=float)
+ 'wetting': np.zeros(self.pressure['wetting'].vector()[:].size, dtype=float),
+ 'nonwetting': np.zeros(self.pressure['nonwetting'].vector()[:].size, dtype=float)
}
for ind in self.has_interface:
interface = self.interface[ind]
@@ -369,7 +346,7 @@ class DomainPatch(df.SubDomain):
neighbour_iter_num = interface.current_iteration[neighbour]
ds = self.ds(ind)
# needed for the read_gli_dofs() functions
- interface_dofs = self._dof_indices_of_interface[ind],#
+ interface_dofs = self._dof_indices_of_interface[ind]
# for each interface, this is a list of dofs indices belonging to another
# interface. The dofs corresponding to these indices must be multiplied
# by 1/2 to to get an average of the gli terms for dofs belonging to
@@ -397,38 +374,36 @@ class DomainPatch(df.SubDomain):
)
if self.isRichards:
- subdomain_has_phases = ['wetting']
# assuming that we normalised atmospheric pressure to zero,
# pc = pw in the Richards case.
pc_prev = pw_prev_iter['wetting']
else:
- subdomain_has_phases = ['wetting', 'nonwetting']
pw_prev = pw_prev_iter['wetting']
pnw_prev = pw_prev_iter['nonwetting']
pc_prev = pnw_prev - pw_prev
- for phase in subdomain_has_phases:
+ for phase in self.has_phases:
# viscosity
- mu_a = mu['phase']
+ mu_a = mu[phase]
# relative permeability
- ka = k['phase']
+ ka = k[phase]
# previous pressure iteration
- pa_prev = pw_prev_iter['phase']
+ pa_prev = pw_prev_iter[phase]
# testfunction
- phi_a = self.testfunction['phase']
+ phi_a = self.testfunction[phase]
if phase == 'wetting':
# 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)
- gl0 = (df.dot(flux, n) - Lambda['phase']*pa_prev)
+ gl0 = (df.dot(flux, n) - Lambda[phase]*pa_prev)
gl0_assembled = df.assemble(dt*gl0*phi_a*ds).get_local()
- # if dofs_in_common_with_other_interfaces['phase'].size == 0
+ # if dofs_in_common_with_other_interfaces[phase].size == 0
# there are no dofs that lie on another interface and we can
# skip the following step.
- if dofs_in_common_with_other_interfaces['phase'].size != 0:
- for common_dof in dofs_in_common_with_other_interfaces['phase']:
+ if dofs_in_common_with_other_interfaces[phase].size != 0:
+ for common_dof in dofs_in_common_with_other_interfaces[phase]:
# from the standpoint of the subdomain we are currently on,
# each dof can belong maximally to two interfaces. On these
# vertices, common to two interfaces, the dofs of the
@@ -440,7 +415,7 @@ class DomainPatch(df.SubDomain):
gl0_assembled[common_dof] = 0.5*gl0_assembled[common_dof]
# now, add the just assembled and corrected term to the Global
# gli_assembled_tmp vector
- gli_assembled_tmp['alpha'] += gl0_assembled
+ gli_assembled_tmp[phase] += gl0_assembled
else:
# phase == 'nonwetting'
# if we are in the two-phase case but our neighbour assumes
@@ -449,23 +424,23 @@ class DomainPatch(df.SubDomain):
# if the neighbour of our subdomain (with index self.subdomain_index)
# assumes a Richards model, we don't have gli terms for the
# nonwetting phase and assemble the terms as zero form.
- gli_assembled_tmp['phase'] += df.assemble(df.Constant(0)*ds).get_local()
+ gli_assembled_tmp[phase] += df.assemble(df.Constant(0)*ds).get_local()
else:
# in this case the neighbour assumes two-phase flow
# and we need to calculte a gl0 torm for the nonwetting
# phase.
# we need anothre flux in this case
flux = -1/mu_a*ka(1-S(pc_prev))*df.grad(pa_prev)
- gl0 = (df.dot(flux, n) - Lambda['phase']*pa_prev)
+ gl0 = (df.dot(flux, n) - Lambda[phase]*pa_prev)
gl0_assembled = df.assemble(dt*gl0*phi_a*ds).get_local()
- if dofs_in_common_with_other_interfaces['phase'].size != 0:
- for common_dof in dofs_in_common_with_other_interfaces['phase']:
+ if dofs_in_common_with_other_interfaces[phase].size != 0:
+ for common_dof in dofs_in_common_with_other_interfaces[phase]:
# see previous case.
gl0_assembled[common_dof] = 0.5*gl0_assembled[common_dof]
# now, add the just assembled and corrected term to the Global
# gli_assembled_tmp vector
- gli_assembled_tmp['phase'] += gl0_assembled
+ gli_assembled_tmp[phase] += gl0_assembled
# writing out glo to the current iteration dictionary of the
# interface takes place after the loop, see below.
@@ -493,12 +468,8 @@ class DomainPatch(df.SubDomain):
interface.gli_term_prev[subdomain].update(#
{'nonwetting': interface.gli_term[neighbour]['nonwetting']}
)
- if self.isRichards:
- subdomain_has_phases = ['wetting']
- else:
- subdomain_has_phases = ['wetting', 'nonwetting']
- for phase in subdomain_has_phases:
+ for phase in self.has_phases:
# in case phase == 'nonwetting' only proceed if the neighbouring
# subdomain does not assume a Richards model.
if (phase == 'wetting') or (not interface.neighbour_isRichards[subdomain]):
@@ -536,13 +507,13 @@ class DomainPatch(df.SubDomain):
gli_prev_neighbour = gli_readout_tmp
gli = gli_p_part + gli_prev_neighbour
# check if there are shared dofs with another interface.
- if dofs_in_common_with_other_interfaces['phase'].size != 0:
- for common_dof in dofs_in_common_with_other_interfaces['phase']:
+ if dofs_in_common_with_other_interfaces[phase].size != 0:
+ for common_dof in dofs_in_common_with_other_interfaces[phase]:
# see previous case.
gli[common_dof] = 0.5*gli[common_dof]
# now, add the just assembled and corrected term to the Global
# gli_assembled_tmp vector
- gli_assembled_tmp['phase'] += gli
+ gli_assembled_tmp[phase] += gli
if neighbour_iter_num -1 == iteration:
# gli_term_prev has been used by the above and can now safely
@@ -570,9 +541,9 @@ class DomainPatch(df.SubDomain):
# save the result of the calculation to the subdomain and to the interface
# dictionaries.
if self.isRichards:
- self.gli_assembled.update({
- 'wetting': gli_assembled_tmp['wetting']
- })
+ self.gli_assembled.update(
+ {'wetting': gli_assembled_tmp['wetting']}
+ )
for ind in self.has_interface:
# self._calc_gli_term_on(interface, iteration)
interface = self.interface[ind]
@@ -595,11 +566,11 @@ class DomainPatch(df.SubDomain):
interface = self.interface[ind]
for phase in ['wetting', 'nonwetting']:
# write gli dofs to interface dicts for communiction
- interface.read_gli_dofs(from_array = self.gli_assembled['phase'], #
- interface_dofs = interface_dofs['phase'],#
- dof_to_vert_map = self.dof2vertex['phase'],#
+ interface.read_gli_dofs(from_array = self.gli_assembled[phase], #
+ interface_dofs = interface_dofs[phase],#
+ dof_to_vert_map = self.dof2vertex[phase],#
local_to_parent_vertex_map = self.parent_mesh_index,#
- phase = 'phase',#
+ phase = phase,#
subdomain_ind = subdomain,#
previous_iter = False
)
@@ -625,15 +596,21 @@ class DomainPatch(df.SubDomain):
"""
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'])
- self.neighbouring_interface_pressure = {'wetting' : df.interpolate(df.Constant(0), self.function_space['wetting'])}
+ self.trialfunction = {
+ 'wetting': df.TrialFunction(self.function_space['wetting'])
+ }
+ self.testfunction = {
+ 'wetting': df.TestFunction(self.function_space['wetting'])
+ }
+ self.neighbouring_interface_pressure = {
+ 'wetting' : df.interpolate(df.Constant(0), self.function_space['wetting'])
+ }
else:
self.function_space = {#
'wetting' : df.FunctionSpace(self.mesh, 'P', 1),#
'nonwetting' : df.FunctionSpace(self.mesh, 'P', 1)#
}
- self.pressure = {
+ self.trialfunction = {
'wetting' : df.TrialFunction(self.function_space['wetting']),
'nonwetting' : df.TrialFunction(self.function_space['nonwetting'])
}
@@ -688,7 +665,7 @@ class DomainPatch(df.SubDomain):
self.interface_marker.set_all(0)
for glob_index in self.has_interface:
# each interface gets marked with the global interface index.
- self.interface[glob_index].mark(self.interface_marker, glob_index)
+ self.interface[glob_index].mark(self.interface_marker, glob_index+1)
# create outerBoundary objects and mark outer boundaries with 1
for index, boundary_points in self.outer_boundary_def_points.items():
# if our domain patch has no outer boundary, this should be reflected
@@ -700,7 +677,7 @@ class DomainPatch(df.SubDomain):
self.outer_boundary.append(#
BoundaryPart(vertices = boundary_points,#
internal = False,
- tol=self.mesh.hmin()/100 )
+ tol= self.tol) #self.mesh.hmin()/100
)
self.outer_boundary[index].mark(self.outer_boundary_marker, 1)
@@ -722,18 +699,26 @@ class DomainPatch(df.SubDomain):
V = self.function_space
marker = self.interface_marker
for ind in self.has_interface:
+ # the marker value for the interfaces is always global interface index+1
+ marker_value = ind + 1
self._dof_indices_of_interface.update(#
{ind: dict()}
)
- if self.isRichards:
- self._dof_indices_of_interface[ind].update(#
- {'wetting': self.interface[ind].dofs_on_interface(V['wetting'], marker, ind)}
- )
- else:
+ for phase in self.has_phases:
self._dof_indices_of_interface[ind].update(#
- {'wetting': self.interface[ind].dofs_on_interface(V['wetting'], marker, ind),#
- 'nonwetting': self.interface[ind].dofs_on_interface(V['nonwetting'], marker, ind)}
+ {phase: self.interface[ind].dofs_on_interface(V[phase], marker, marker_value)}
)
+ print(f"subdomain{self.subdomain_index}: dofs on interface {ind}:", self._dof_indices_of_interface[ind][phase])
+ # if self.isRichards:
+ # self._dof_indices_of_interface[ind].update(#
+ # {'wetting': self.interface[ind].dofs_on_interface(V['wetting'], marker, marker_value)}
+ # )
+ # print(f"subdomain{self.subdomain_index}: dofs on interface {ind}:", self._dof_indices_of_interface[ind]['wetting'])
+ # else:
+ # self._dof_indices_of_interface[ind].update(#
+ # {'wetting': self.interface[ind].dofs_on_interface(V['wetting'], marker, marker_value),#
+ # 'nonwetting': self.interface[ind].dofs_on_interface(V['nonwetting'], marker, marker_value)}
+ # )
def _calc_interface_has_common_dof_indices(self):
""" populate dictionary self._interface_has_common_dof_indices
@@ -944,10 +929,10 @@ class BoundaryPart(df.SubDomain):
# the test ((p[0] - xmax)*(p[0] - xmin) < 0) might fail if p[0] and one of p1[0]
# or p2[0] are equal. In pp1 or pp2 is vertical. We need to still calculate
# the distance to the line segment in this case
- if ((p[0] - xmax)*(p[0] - xmin) < 0) or np.absolute((p[0] - xmax)*(p[0] - xmin)) < tol:
+ if ((p[0] - xmax)*(p[0] - xmin) < tol) or np.absolute((p[0] - xmax)*(p[0] - xmin)) < tol:
# here we have a chance to actually lie on the line segment.
segment_direction = (p2 - p1)/np.linalg.norm(p2 - p1)
- distance = (p - p1) - np.dot(segment_direction, p - p1) * segment_direction
+ distance = (p - p1) - np.dot(segment_direction, p - p1)*segment_direction
# check again for equality with p1 and p2 just to be sure. =)
if (np.linalg.norm(p - p1) < tol) or (np.linalg.norm(p - p2) < tol):
# print(f"grid node at {p} was close to either {p1} or {p2}")
@@ -1069,7 +1054,9 @@ class Interface(BoundaryPart):
according to self.inside
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,#
+ print_vertex_indices = False)
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
diff --git a/RR-2-patch-test-case/RR-2-patch-test.py b/RR-2-patch-test-case/RR-2-patch-test.py
index 55af1aace4ce587fb8d20fb51cee5031b4ad52f0..fd6c65c885253c0e93126330677a2e70c85e648a 100755
--- a/RR-2-patch-test-case/RR-2-patch-test.py
+++ b/RR-2-patch-test-case/RR-2-patch-test.py
@@ -183,7 +183,7 @@ dirichletBC = exact_solution
mesh_resolution = 3
# initialise LDD simulation class
-simulation = ldd.LDDsimulation()
+simulation = ldd.LDDsimulation(tol = 1E-12)
simulation.set_parameters(output_dir = "",#
subdomain_def_points = subdomain_def_points,#
isRichards = isRichards,#
@@ -211,7 +211,7 @@ domain_marker = simulation.domain_marker
mesh_subdomain = simulation.mesh_subdomain
# mesh = mesh_subdomain[0]
-# interface_marker = df.MeshFunction('size_t', mesh, mesh.topology().dim()-1)
+# interface_marker = df.MeshFunction('int', mesh, mesh.topology().dim()-1)
# interface_marker.set_all(0)
# interface = dp.Interface(vertices=interface_def_points[0], #
# tol=mesh.hmin()/100,#
@@ -226,53 +226,68 @@ interface_marker = simulation.interface_marker
subdoms = simulation.subdomain
-for iterations in range(1):
- #
-
- a1 = (L1*u1*v1)*dx1 + (df.inner(relative_permeability(saturation(p1_i, 1), 1)*df.grad(u1), df.grad(v1)))*dx_1 + (timestep*lambda1*u1*v1)*ds1(1) #timestep*
- rhs1 = (L1*p1_i*v1)*dx1 - ((saturation(p1_i, 1) - saturation(p1_0, 1))*v1)*dx1 + (timestep*(source1 - g1_i)*v1)*ds1(1)
- rhssplit1 = (L1*p1_i*v1)*dx1 - ((saturation(p1_i, 1) - saturation(p1_0, 1))*v1)*dx1
- extratrem = (timestep*(source1 - g1_i)*v1)*ds1(1)
- splitrhs = df.assemble(rhssplit1)
- print("splitrhs: \n", splitrhs.get_local())
- extrat_assembled = df.assemble(extratrem)
- print("extratrem: \n", extrat_assembled.get_local())
- added = splitrhs + extrat_assembled
- print("both added together:\n", added.get_local())
-
- A1 = df.assemble(a1)
- dsform = (timestep*lambda1*u1*v1)*ds1(1)
- dsform_vec = df.assemble(dsform)
- b1 = df.assemble(rhs1)
- #p_gamma1.apply(A1,b1)
- # outerBC1.apply(A1,b1)
- u = df.Function(V1)
- U1 = u.vector()
- df.solve(A1,U1,b1)
- # print("Form:\n", A1.array())
- print("RHS:\n", b1.get_local())
- #print("ds term:\n", dsform_vec.array())
- # print('solution:\n', U1.get_local())
- p1_i.vector()[:] = U1
- interface[0].read_pressure_dofs(from_function = p1_i, #
- interface_dofs = dofs_on_interface1,#
- dof_to_vert_map = dom1_d2v,#
- local_to_parent_vertex_map = dom1_parent_vertex_indices,#
- phase = 'wetting',#
- subdomain_ind = 1)
-#
-# Save mesh to file
-df.File('./test_domain_layered_soil.xml.gz') << mesh_subdomain[0]
-df.File('./test_domain_mesh.pvd') << mesh_subdomain[0]
-df.File('./test_global_interface_marker.pvd') << interface_marker
-df.File('./test_subdomain1.xml.gz') << mesh_subdomain[1]
-df.File('./test_subdomain2.xml.gz') << mesh_subdomain[2]
-df.File('./test_domain_marker.pvd') << domain_marker
-df.File('./test_domain_layered_soil_solution.pvd') << u
-df.File('./test_subdomain1_interface_marker.pvd') << interface_marker1
-df.File('./test_subdomain2_interface_marker.pvd') << interface_marker2
-df.File('./test_subdomain1_boundary_marker.pvd') << boundary_marker1
-df.File('./test_subdomain2_boundary_marker.pvd') << boundary_marker2
+# df.File('./test_domain_layered_soil.xml.gz') << mesh_subdomain[0]
+# df.File('./test_domain_mesh.pvd') << mesh_subdomain[0]
+df.File('./global_interface_marker.pvd') << interface_marker
+# df.File('./test_subdomain1.xml.gz') << mesh_subdomain[1]
+# df.File('./test_subdomain2.xml.gz') << mesh_subdomain[2]
+df.File('./domain_marker.pvd') << domain_marker
+# df.File('./test_domain_layered_soil_solution.pvd') << u
+df.File('./subdomain1_interface_marker.pvd') << simulation.subdomain[1].interface_marker
+df.File('./subdomain2_interface_marker.pvd') << simulation.subdomain[2].interface_marker
+# df.File('./test_subdomain1_boundary_marker.pvd') << boundary_marker1
+# df.File('./test_subdomain2_boundary_marker.pvd') << boundary_marker2
+
+simulation.LDDsolver(time = 0)
+# df.info(parameters, True)
+
+# for iterations in range(1):
+# #
+#
+# a1 = (L1*u1*v1)*dx1 + (df.inner(relative_permeability(saturation(p1_i, 1), 1)*df.grad(u1), df.grad(v1)))*dx_1 + (timestep*lambda1*u1*v1)*ds1(1) #timestep*
+# rhs1 = (L1*p1_i*v1)*dx1 - ((saturation(p1_i, 1) - saturation(p1_0, 1))*v1)*dx1 + (timestep*(source1 - g1_i)*v1)*ds1(1)
+# rhssplit1 = (L1*p1_i*v1)*dx1 - ((saturation(p1_i, 1) - saturation(p1_0, 1))*v1)*dx1
+# extratrem = (timestep*(source1 - g1_i)*v1)*ds1(1)
+# splitrhs = df.assemble(rhssplit1)
+# print("splitrhs: \n", splitrhs.get_local())
+# extrat_assembled = df.assemble(extratrem)
+# print("extratrem: \n", extrat_assembled.get_local())
+# added = splitrhs + extrat_assembled
+# print("both added together:\n", added.get_local())
+#
+# A1 = df.assemble(a1)
+# dsform = (timestep*lambda1*u1*v1)*ds1(1)
+# dsform_vec = df.assemble(dsform)
+# b1 = df.assemble(rhs1)
+# #p_gamma1.apply(A1,b1)
+# # outerBC1.apply(A1,b1)
+# u = df.Function(V1)
+# U1 = u.vector()
+# df.solve(A1,U1,b1)
+# # print("Form:\n", A1.array())
+# print("RHS:\n", b1.get_local())
+# #print("ds term:\n", dsform_vec.array())
+# # print('solution:\n', U1.get_local())
+# p1_i.vector()[:] = U1
+# interface[0].read_pressure_dofs(from_function = p1_i, #
+# interface_dofs = dofs_on_interface1,#
+# dof_to_vert_map = dom1_d2v,#
+# local_to_parent_vertex_map = dom1_parent_vertex_indices,#
+# phase = 'wetting',#
+# subdomain_ind = 1)
+# #
+# # Save mesh to file
+# df.File('./test_domain_layered_soil.xml.gz') << mesh_subdomain[0]
+# df.File('./test_domain_mesh.pvd') << mesh_subdomain[0]
+# df.File('./test_global_interface_marker.pvd') << interface_marker
+# df.File('./test_subdomain1.xml.gz') << mesh_subdomain[1]
+# df.File('./test_subdomain2.xml.gz') << mesh_subdomain[2]
+# df.File('./test_domain_marker.pvd') << domain_marker
+# df.File('./test_domain_layered_soil_solution.pvd') << u
+# df.File('./test_subdomain1_interface_marker.pvd') << interface_marker1
+# df.File('./test_subdomain2_interface_marker.pvd') << interface_marker2
+# df.File('./test_subdomain1_boundary_marker.pvd') << boundary_marker1
+# df.File('./test_subdomain2_boundary_marker.pvd') << boundary_marker2
# boundary_marker1 = simulation.subdomain[1].outer_boundary_marker
# boundary_marker2 = simulation.subdomain[2].outer_boundary_marker
@@ -281,8 +296,8 @@ df.File('./test_subdomain2_boundary_marker.pvd') << boundary_marker2
# ds_1 = simulation.subdomain[1].ds
# ds_2 = simulation.subdomain[2].ds
#
-# interface_marker1 = simulation.subdomain[1].interface_marker#df.MeshFunction('size_t', mesh_subdomain[1], mesh_subdomain[1].topology().dim()-1)
-# interface_marker2 = simulation.subdomain[2].interface_marker#df.MeshFunction('size_t', mesh_subdomain[2], mesh_subdomain[2].topology().dim()-1)
+# interface_marker1 = simulation.subdomain[1].interface_marker#df.MeshFunction('int', mesh_subdomain[1], mesh_subdomain[1].topology().dim()-1)
+# interface_marker2 = simulation.subdomain[2].interface_marker#df.MeshFunction('int', mesh_subdomain[2], mesh_subdomain[2].topology().dim()-1)
# # interface_marker1.set_all(0)
# # interface_marker2.set_all(0)
# # print(dir(interface_marker1))