diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index 87fbf67557192ac4b43bb93760cb9b0ad78a9d9d..84c441f4266ec2c377704e8be09e8d55ef109c8d 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -588,32 +588,19 @@ class LDDsimulation(object):
"""
subdomain = self.subdomain[subdomain_index]
iteration = subdomain.iteration_number
- # 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 = phase)
form = governing_problem['form']
- rhs_without_gli = governing_problem['rhs_without_gli']
+ rhs = governing_problem['rhs']
# assemble the form and rhs
form_assembled = df.assemble(form)
- rhs_without_gli_assembled = df.assemble(rhs_without_gli)
- # 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.
- 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())
- # we need the rhs as a dolfin function, so we overwrite
- # rhs_without_gli_assembled with the gli_term subtracted from itself.
- rhs_without_gli_assembled.set_local(rhs_without_gli_assembled.get_local() - gli_term_assembled)
- rhs_assembled = rhs_without_gli_assembled
+ rhs_assembled = df.assemble(rhs)
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}: rhs:\n", rhs_assembled.get_local())
print(f"phase = {phase}: gli_term:\n", gli_term_assembled)
print(f"phase = {phase}: rhs_assembled:\n", rhs_assembled.get_local())
@@ -627,12 +614,6 @@ class LDDsimulation(object):
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())
linear_solver = subdomain.linear_solver
@@ -644,6 +625,7 @@ class LDDsimulation(object):
if debug and subdomain.mesh.num_cells() < 36:
print("\npressure after solver:\n", subdomain.pressure[phase].vector().get_local())
+
def _init_solution_files(self):
""" set up solution files for saving the solution of the LDD scheme for
each subdomain.
diff --git a/LDDsimulation/domainPatch.py b/LDDsimulation/domainPatch.py
index 989b79f555a2d81f9f11107162bd75d3b483b68c..1d01e0d49f068af7f55cc49bbbc9c64b784754dd 100644
--- a/LDDsimulation/domainPatch.py
+++ b/LDDsimulation/domainPatch.py
@@ -339,22 +339,26 @@ class DomainPatch(df.SubDomain):
if phase == "wetting":
# gather interface forms
interface_forms = []
+ gli_forms = []
# the marker values of the interfacese are global interfaces index + 1
for interface in self.has_interface:
# print(f"subdomain{self.subdomain_index} has interfaces: {interface}\n")
marker_value = self.interface[interface].marker_value
# interface_forms += (dt*Lambda_a*pa*phi_a)*ds(interface)
interface_forms.append((dt*Lambda_a*pa*phi_a)*ds(marker_value))
+ gli = self.calc_gli_term(interface_index=interface, phase=phase)
+ gli_forms.append((dt*gli*phi_a)*ds(marker_value))
# 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
+ 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 = []
+ gli_forms = []
for interface in self.has_interface:
marker_value = self.interface[interface].marker_value
if self.interface[interface].neighbour_isRichards[self.subdomain_index]:
@@ -364,15 +368,18 @@ class DomainPatch(df.SubDomain):
# interface_forms += (df.Constant(0)*phi_a)*ds(interface)
# pass
interface_forms.append((0*pa*phi_a)*ds(marker_value))
+ gli_forms.append((0*phi_a)*ds(marker_value))
else:
# interface_forms += (dt*Lambda_a*pa*phi_a)*ds(interface)
interface_forms.append((dt*Lambda_a*pa*phi_a)*ds(marker_value))
+ gli = self.calc_gli_term(interface_index=interface, phase=phase)
+ gli_forms.append((dt*gli*phi_a)*ds(marker_value))
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
+ 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 ')
@@ -380,13 +387,13 @@ class DomainPatch(df.SubDomain):
form = form1 + form2 + sum(interface_forms)
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
+ rhs = rhs1 + rhs2 + dt*source_a*phi_a*dx - sum(gli_forms) - rhs_gravity
else:
- rhs_without_gli = rhs1 + rhs2 + dt*source_a*phi_a*dx
+ rhs = rhs1 + rhs2 + dt*source_a*phi_a*dx - sum(gli_forms)
form_and_rhs = {#
'form': form,#
- 'rhs_without_gli' : rhs_without_gli
+ 'rhs' : rhs
}
return form_and_rhs
@@ -394,7 +401,8 @@ class DomainPatch(df.SubDomain):
def calc_gl0_term(self, debug: bool = False) -> None: # tp.Dict[str, fl.Form]
"""calculate the gl0 terms for all interfaces of the subdomain.
- assemble the gl0 terms for all interfaces of the subdomain and save them
+ calculate the gl0 terms for all interfaces of the subdomain, interpolate
+ them to our functionspace and save them
to the dictionries of the interfaces. This method gets called by the LDDsolver
before entering the L-iterations in each time step.
This method assumes that self.pressure_prev_timestep has assigned the values
@@ -403,36 +411,9 @@ class DomainPatch(df.SubDomain):
subdomain = self.subdomain_index
dt = self.timestep_size
Lambda = self.lambda_param
- # since we loop over all interfaces in self.has_interface we need a temporary
- # vectors (numpy arrays) for each phase that we will use to add all dofs of the gli terms
- # into one array.
- gli_assembled_tmp = dict()
- for phase in self.has_phases:
- # print("number of dofs: ", self.pressure_prev_iter['wetting'].vector()[:].size)
- gli_assembled_tmp.update(
- {phase: np.zeros(self.pressure[phase].vector()[:].size, dtype=float)}
- )
for interf_ind in self.has_interface:
interface = self.interface[interf_ind]
- # # neighbour index
- # neighbour = interface.neighbour[subdomain]
- # neighbour_iter_num = interface.current_iteration[neighbour]
- ds = self.ds(interface.marker_value)
- # needed for the read_gli_dofs() functions
- interface_dofs = self._dof_indices_of_interface[interf_ind]
- # for each interface, the following 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
- # two different interfaces.
- dofs_in_common_with_other_interfaces = self._interface_has_common_dof_indices[interf_ind]
- if debug:
- neighbour = interface.neighbour[subdomain]
- print(f"On subdomain{subdomain}, we have:\n",
- f"Interface{interf_ind}, Neighbour index = {neighbour}\n",
- f"dofs on interface:\n{interface_dofs['wetting']}\n",
- #f"dof2global_vertex_map:\n{self.parent_mesh_index[self.dof2vertex['wetting'][interface_dofs['wetting']]]}\n"
- f"and dofs common with other interfaces:\n{dofs_in_common_with_other_interfaces['wetting']}")
if self.isRichards:
# assuming that we normalised atmospheric pressure to zero,
@@ -453,7 +434,7 @@ class DomainPatch(df.SubDomain):
# previous timestep pressure
pa_prev = self.pressure_prev_timestep[phase]
# testfunction
- phi_a = self.testfunction[phase]
+ # phi_a = self.testfunction[phase]
if self.include_gravity:
z_a = self.gravity_term[phase]
if phase == 'nonwetting' and interface.neighbour_isRichards[subdomain]:
@@ -486,257 +467,141 @@ class DomainPatch(df.SubDomain):
# part of gl0 seperately at the intersection points of two
# interfaces (more details in comment below). Therefor we save
# the robin pressure part and the neumann flux part separately
- robin_pressure = Lambda[phase]*pa_prev
- robin_pressure_assembled = df.assemble(dt*robin_pressure*phi_a*ds).get_local()
- neumann_flux_assembled = df.assemble(dt*df.dot(flux, n)*phi_a*ds).get_local()
- # gl0_assembled = df.assemble(dt*gl0*phi_a*ds).get_local()
-
- # 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 we have dofs in common with another interface we get two
- # different gl0 terms from both interfaces apriori, so in
- # order to get one single dof for the intersection we do:
- # we replace the neuman flux by the flux in dircetion of the
- # linear combination of both adjacent normal vectors:
- # flux_neumann = df.dot(flux, n_1 + 1_2)/norm(n_1 + 1_2)
- # since we are looping over the interfaces we correct each dof
- # that is common with another interfaces with the corresponding
- # factor norm(n_1 + 1_2).
- # We take the mean of the pressure on both interfaces, so
- # we multiply these by 1/2
- if dofs_in_common_with_other_interfaces[phase].size != 0:
- if debug:
- print(f"\nthe following dofs (phase = {phase}) of interface{interf_ind} are in common with \n",
- f" other interfaces: {dofs_in_common_with_other_interfaces[phase]}")
- for common_dof in dofs_in_common_with_other_interfaces[phase]:
- n_norm = self._normal_vector_norms_for_common_dofs[interf_ind][phase][common_dof]
- # from the standpoint of the subdomain we are currently on,
- # each dof can belong maximally to two interfaces.
- print(f"On subdomain{self.subdomain_index} and interface{interf_ind} the flux\n",
- f"neumann_flux_assembled[{common_dof}] is devided by the norm {n_norm}",
- f"coordinate of the common_dof = {self.mesh.coordinates()[self.dof2vertex[phase][common_dof]]}")
- neumann_flux_assembled[common_dof] = neumann_flux_assembled[common_dof]/n_norm
- robin_pressure_assembled[common_dof] = robin_pressure_assembled[common_dof]/2
-
- gl0_assembled = neumann_flux_assembled - robin_pressure_assembled
- # now, add the just assembled and corrected term to the Global
- # gli_assembled_tmp vector
- gli_assembled_tmp[phase] += gl0_assembled
-
- # writing out glo to the current iteration dictionary of the
- # interface takes place after the loop, see below.
-
- ### END for ind in self.has_interface:
- # save the result of the calculation to the subdomain and to the interface
- # dictionaries.
- for phase in self.has_phases:
- for interf_ind in self.has_interface:
- interface = self.interface[interf_ind]
- # # needed for the read_gli_dofs() functions
- interface_dofs = self._dof_indices_of_interface[interf_ind]
- # if debug:
- # print(f"Interface{interf_ind}.gli_term_prev[{subdomain}][{phase}]=\n",interface.gli_term_prev[subdomain][phase])
- # print(f"before writing to interface[{interf_ind}] the values gli_assembled[{phase}]=\n",gli_assembled_tmp[phase])
- # write gli dofs to interface dicts for communiction.
- interface.read_gli_dofs(from_array = gli_assembled_tmp[phase],#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,#
- subdomain_ind = subdomain,#
- previous_iter = True
- )
- # if debug:
- # print(f"gl0_terms after writing to interface[{interf_ind}] Interface[{interf_ind}].gli_term_prev[{subdomain}][{phase}]=\n",interface.gli_term_prev[subdomain][phase])
- ### END calc_gl0_term
+ robin_pressure = Lambda[phase]*pa_prev.vector()
+ # normal_flux = df.dot(flux, n).get_local()
+ # print("type(normal_flux)", type(normal_flux))
+ V = self.function_space[phase]
+ degree = V.ufl_element().degree()
+ interface_mesh = df.restriction(self.interface_marker, interface.marker_value)
+ # W = VectorFunctionSpace(interface_mesh, 'Q', 0)
+ W = FunctionSpace(interface_mesh, 'Q', 0)
+ # Project the normal flux onto V so that we can
+ normal_flux_Q0 = df.project(df.dot(flux,n), W)
+ flux_x, flux_y = flux_FEM.split(deepcopy=True)
+ normal_flux = normal_flux_FEM.vector()
+ interface_dofs = self._dof_indices_of_interface[interf_ind][phase]
+ parent = self.parent_mesh_index
+ d2v = self.dof2vertex[phase]
+ for dof_index in interface_dofs:
+ gli = normal_flux[dof_index] - robin_pressure[dof_index]
+ # from_array 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.
+ interface.gli_term_prev[subdomain_ind][phase].update(
+ {parent[d2v[dof_index]]: gli}
+ )
+ ### END calc_gl0_term
- def calc_gli_term(self, debug: bool = False) -> None: # tp.Dict[str, fl.Form]
- """calculate the gli terms for the iteration number of the subdomain
- calculate the gl terms for the iteration number of the subdomain and save them
- to self.gli_assembled. The term gets saved as an already assembled sparse
- ds form to be added to the rhs by the solver. This is due to the need of
- having to communicate dofs over the interface. Additionally, the dofs are
- being saved to an interface dictionary exactly as the pressures.
+ def calc_gli_term(self,
+ interface_index: int,
+ phase: str,
+ debug: bool = False) -> df.Function: # tp.Dict[str, fl.Form]
+ """calculate the gli term for phase phase, interface interface_index
+ and the iteration number of the subdomain
- Caution: The array self.gli_assembled needs to be added to the rhs with
- a minus sign by the solver!
+ calculate the gl term for phase phase, interface interface_index and
+ the iteration number of the subdomain and return it as a dolfin Function
+ used for assembling the governing form in self.govering_problem. Additionally,
+ save the newly computed gli_term to the gli dictionries of the interface.
"""
# this is the number of the current iteration of the subdomain.
iteration = self.iteration_number
subdomain = self.subdomain_index
- dt = self.timestep_size
- Lambda = self.lambda_param
- # since we loop over all interfaces in self.has_interface we need a temporary
- # vector (numpy array) that we will use to add all dofs of the gli terms
- # into one array.
- gli_assembled_tmp = dict()
- for phase in self.has_phases:
- # print("number of dofs: ", self.pressure_prev_iter['wetting'].vector()[:].size)
- gli_assembled_tmp.update(
- {phase: np.zeros(self.pressure[phase].vector()[:].size, dtype=float)}
+ interface = self.interface[interface_index]
+ # gli should be initialised by zero.
+ gli = df.Function(self.function_space[phase])
+ # neighbour index
+ neighbour = interface.neighbour[subdomain]
+ # update the current_iteration number of subdomain stored in the interface
+ # to the current iteration number of the subdomain.
+ interface.current_iteration[subdomain] = iteration
+ # save the iteration number of the neighbour
+ neighbour_iter_num = interface.current_iteration[neighbour]
+ # needed for the read_gli_dofs() functions
+ interface_dofs = self._dof_indices_of_interface[interface_index][phase]
+ # for each interface, the following 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
+ # two different interfaces.
+ dofs_in_common_with_other_interfaces = self._interface_has_common_dof_indices[interface_index]
+ if debug:
+ print(f"On subdomain{subdomain}, we have:\n",
+ f"Interface{interface_index}, Neighbour index = {neighbour}\n",
+ f"dofs on interface:\n{interface_dofs['wetting']}\n",
+ f"dof2global_vertex_map:\n{self.parent_mesh_index[self.dof2vertex['wetting'][interface_dofs['wetting']]]}\n"
+ f"and dofs common with other interfaces:\n{dofs_in_common_with_other_interfaces['wetting']}")
+
+ if (phase == 'wetting') or (not interface.neighbour_isRichards[subdomain]):
+ # read previous neighbouring pressure values from interface.
+ # depending on our iteration number and the iteration number of
+ # the neighbour, weed to either read from the previous pressure
+ # dictionary or the current one.
+ if iteration == neighbour_iter_num + 1:
+ # in this case the neighbour has not yet iterated beyond
+ # the iteration we are currently in, so
+ # interface.pressure_values[neighbour] holds the "previous pressure"
+ # we need to read from the interface to calculate the current gli.
+ pressure_prev = interface.pressure_values[neighbour][phase]
+
+ # choose dictionary to save the newly calculated dofs to.
+ # in this case, the neighbour has not done the next iteration
+ # and will need both pressure and gli_term_prev of the current
+ # subdomain. Therefor, we save the the calculated gli_terms to
+ # the gli_term dictionary of the current subdomain, not the
+ # gli_term_prev dictionary.
+ gli_save_dictionary = interface.gli_term[subdomain][phase]
+
+ elif iteration == neighbour_iter_num:
+ # this is the only possible other case. In this case the
+ # neighbour has iterated once beyond the iteration we are
+ # currently in, so interface.pressure_values_prev holds the "previous pressure"
+ # we need to read from the interface to calculate gli.
+ pressure_prev = interface.pressure_values_prev[neighbour][phase]
+
+ # choose dictionary to save the newly calculated dofs to.
+ # This should be the case iteration == neighbour_iter_num:
+ # in this case, the neighbour has done the next iteration already
+ # and will not need both pressure_prev and gli_term_prev of the current
+ # subdomain. Therefor, we save the the calculated gli_terms to
+ # the gli_term_prev dictionary of the current subdomain, not the
+ # gli_term dictionary.
+ gli_save_dictionary = interface.gli_term_prev[subdomain][phase]
+ else:
+ raise RuntimeError("\n!!!!!!!!! Warning !!!!!!!!!!!!!! \n"+\
+ "this case should not appear\n"+"neighbour_iter_num =" + "{}".format(neighbour_iter_num)\
+ + "for neighbour = "+ "{}".format(neighbour)+\
+ "and iteration = " + "{} on subdomain".format(iteration)+ "{}".format(subdomain)\
+ + "You suck! Revisite DomainPatch.calc_gli_term() and LDDsimulation.LDDsolver ASAP.")
+
+ parent = self.parent_mesh_index
+ d2v = self.dof2vertex[phase]
+ # read gli_prev term from the neighbour.
+ gli_prev = interface.gli_term_prev[neighbour][phase]
+ dt = self.timestep_size
+ Lambda = self.lambda_param[phase]
+ for dof_index in interface_dofs:
+ glob_vert_ind = parent[d2v[dof_index]]
+ gli_value = -2*Lambda*pressure_prev[glob_vert_ind] - gli_prev[glob_vert_ind]
+ gli.vector()[dof_index] = gli_value
+ gli_save_dictionary[glob_vert_ind] = gli_value
+
+
+ # In the following case, the gli_term of the neighbour is saved to
+ # gli_term currently and will be overwritten in the next step,
+ # if we don't save it to gli_term_prev of the neighbour.
+ # similarly for the pressure values.
+ if iteration == neighbour_iter_num:
+ # print(f"we are on subdomain{subdomain} and interface{interface.global_index} (alt index:{ind}) has neighbouring \n",
+ # f"subdomains: subdomain{subdomain} and subdomain{neighbour}")
+ interface.gli_term_prev[neighbour].update(#
+ {phase: interface.gli_term[neighbour][phase].copy()}
+ )
+ interface.pressure_values_prev[neighbour][phase].update(
+ {phase: interface.pressure_values[neighbour][phase].copy()}
)
- for ind in self.has_interface:
- interface = self.interface[ind]
- # neighbour index
- neighbour = interface.neighbour[subdomain]
- # update the current_iteration number of subdomain stored in the interface
- # to the current iteration number of the subdomain.
- interface.current_iteration[subdomain] = iteration
- # save the iteration number of the neighbour
- neighbour_iter_num = interface.current_iteration[neighbour]
- ds = self.ds(interface.marker_value)
- # needed for the read_gli_dofs() functions
- interface_dofs = self._dof_indices_of_interface[ind]
- # for each interface, the following 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
- # two different interfaces.
- dofs_in_common_with_other_interfaces = self._interface_has_common_dof_indices[ind]
- if debug:
- print(f"On subdomain{subdomain}, we have:\n",
- f"Interface{ind}, Neighbour index = {neighbour}\n",
- f"dofs on interface:\n{interface_dofs['wetting']}\n",
- f"dof2global_vertex_map:\n{self.parent_mesh_index[self.dof2vertex['wetting'][interface_dofs['wetting']]]}\n"
- f"and dofs common with other interfaces:\n{dofs_in_common_with_other_interfaces['wetting']}")
-
- 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]):
- # read gli_prev term from the neighbour.
- # to be able to sum the gli terms of all interfaces together,
- # we need a dummy vector to hold the dofs read from the interface.
- gli_prev_of_neighbour = np.zeros(gli_assembled_tmp[phase].size, dtype = float)
- interface.write_gli_dofs(to_array = gli_prev_of_neighbour, #
- interface_dofs = interface_dofs[phase],#
- dof_to_vert_map = self.dof2vertex[phase],#
- local_to_parent_vertex_map = self.parent_mesh_index,#
- phase = phase,#
- subdomain_ind = neighbour,#
- # this is set to True because we need gli_prev
- previous_iter = True
- )
- if debug:
- print(f"read gli_prev_iter[{subdomain}]:\n", gli_prev_of_neighbour)
-
- # read previous neighbouring pressure values from interface and write them to
- # the function pa_prev. We need to make a distinction first.
- if iteration == neighbour_iter_num + 1:
- # in this case the neighbour has not yet iterated beyond
- # the iteration we are currently in, so
- # interface.pressure_values[neighbour] holds the "previous pressure"
- # we need to read from the interface to calculate the current gli.
- take_previous_pressure_from_interface = False
-
- elif iteration == neighbour_iter_num:
- # this is the only possible other case. In this case the
- # neighbour has iterated once beyond the iteration we are
- # currently in, so interface.pressure_values_prev holds the "previous pressure"
- # we need to read from the interface to calculate gli.
- take_previous_pressure_from_interface = True
- else:
- raise RuntimeError("\n!!!!!!!!! Warning !!!!!!!!!!!!!! \n"+\
- "this case should not appear\n"+"neighbour_iter_num =" + "{}".format(neighbour_iter_num)\
- + "for neighbour = "+ "{}".format(neighbour)+\
- "and iteration = " + "{} on subdomain".format(iteration)+ "{}".format(subdomain)\
- + "You suck! Revisite DomainPatch.calc_gli_term() and LDDsimulation.LDDsolver ASAP.")
-
- pa_prev = df.interpolate(df.Constant(0), self.function_space[phase])
- interface.write_pressure_dofs(to_function = pa_prev,#self.neighbouring_interface_pressure[phase], #
- interface_dofs = interface_dofs[phase],#
- dof_to_vert_map = self.dof2vertex[phase],#
- local_to_parent_vertex_map = self.parent_mesh_index,#
- phase = phase,#
- subdomain_ind = neighbour,#
- previous_iter = take_previous_pressure_from_interface)
-
- # use the above to calculate the gli update with
- # in the paper p^{i-1}_{3-l}
- # pa_prev = self.neighbouring_interface_pressure[phase]
- phi_a = self.testfunction[phase]
- gli_pressure_part = df.assemble(-2*dt*Lambda[phase]*pa_prev*phi_a*ds)
- gli_pressure_part_array = gli_pressure_part.get_local()
- # 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]:
- # n_norm = self._normal_vector_norms_for_common_dofs[ind][phase][common_dof]
- gli_prev_of_neighbour[common_dof] = gli_prev_of_neighbour[common_dof]/2
-
- # see previous case.
- gli_pressure_part_array[common_dof] = 0.5*gli_pressure_part_array[common_dof]
-
- gli = gli_pressure_part_array - gli_prev_of_neighbour
- # now, add the just assembled and corrected term to the Global
- # gli_assembled_tmp vector
- gli_assembled_tmp[phase] += gli
-
- ### END for ind in self.has_interface:
- # save the result of the calculation to the subdomain and to the interface
- # dictionaries.
- for phase in self.has_phases:
- self.gli_assembled.update(
- {phase: gli_assembled_tmp[phase].copy()}
- )
- # communicate the newly calculated terms
- for ind in self.has_interface:
- # self._calc_gli_term_on(interface, iteration)
- interface = self.interface[ind]
- neighbour = interface.neighbour[subdomain]
- neighbour_iter_num = interface.current_iteration[neighbour]
- # needed for the read_gli_dofs() functions
- interface_dofs = self._dof_indices_of_interface[ind]
-
- if debug:
- print(f"Interface{ind}.gli_term[{subdomain}][{phase}]=\n",interface.gli_term[subdomain][phase])
- print(f"before writing to interface[{ind}] gli_assembled[{phase}]=\n",self.gli_assembled[phase])
- # write gli dofs to interface dicts for communiction. We have to
- # distinguish different cases.
- if iteration == neighbour_iter_num + 1:
- # in this case, the neighbour has not done the next iteration
- # and will need both pressure and gli_term_prev of the current
- # subdomain. Therefor, we save the the calculated gli_terms to
- # the gli_term dictionary of the current subdomain, not the
- # gli_term_prev dictionary.
- 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,#
- subdomain_ind = subdomain,#
- previous_iter = False
- )
- else:
- # This should be the case iteration == neighbour_iter_num:
- # in this case, the neighbour has done the next iteration already
- # and will not need both pressure_prev and gli_term_prev of the current
- # subdomain. Therefor, we save the the calculated gli_terms to
- # the gli_term_prev dictionary of the current subdomain, not the
- # gli_term dictionary.
- 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,#
- subdomain_ind = subdomain,#
- previous_iter = True
- )
- # Additionally, the gli_term of the neighbour is saved to
- # gli_term currently and will be overwritten in the next step,
- # if we don't save it to gli_term_prev of the neighbour.
- # similarly for the pressure values.
- if phase == 'wetting' or not interface.neighbour_isRichards[subdomain]:
- # print(f"we are on subdomain{subdomain} and interface{interface.global_index} (alt index:{ind}) has neighbouring \n",
- # f"subdomains: subdomain{subdomain} and subdomain{neighbour}")
- interface.gli_term_prev[neighbour].update(#
- {phase: interface.gli_term[neighbour][phase].copy()}
- )
- interface.pressure_values_prev[neighbour][phase].update(
- {phase: interface.pressure_values[neighbour][phase].copy()}
- )
+ return gli
### END calc_gli_term
def null_all_interface_iteration_numbers(self) -> None:
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 f3d946e7b7de6046c95548681e567ca5727eea92..14c9635e712a0cb965a45e59f4bf403552ad55a3 100755
--- a/RR-2-patch-test-case/RR-2-patch-test.py
+++ b/RR-2-patch-test-case/RR-2-patch-test.py
@@ -11,41 +11,42 @@ import functools as ft
##### Domain and Interface ####
# global simulation domain domain
-sub_domain0_vertices = [df.Point(0.0,-1.0), #
+sub_domain0_vertices = [df.Point(-1.0,-1.0), #
df.Point(1.0,-1.0),#
df.Point(1.0,1.0),#
- df.Point(0.0,1.0)]
+ df.Point(-1.0,1.0)]
# interface between subdomain1 and subdomain2
-interface12_vertices = [df.Point(0.0, 0.0),
+interface12_vertices = [df.Point(-1.0, 0.0),
df.Point(1.0, 0.0) ]
# subdomain1.
sub_domain1_vertices = [interface12_vertices[0],
interface12_vertices[1],
- df.Point(1.0,1.0),
- df.Point(0.0,1.0) ]
+ sub_domain0_vertices[2],
+ sub_domain0_vertices[3] ]
# vertex coordinates of the outer boundaries. If it can not be specified as a
# polygon, use an entry per boundary polygon. This information is used for defining
# the Dirichlet boundary conditions. If a domain is completely internal, the
# dictionary entry should be 0: None
subdomain1_outer_boundary_verts = {
- 0: [interface12_vertices[1], #
- df.Point(1.0,1.0), #
- df.Point(0.0,1.0), #
+ 0: [interface12_vertices[1],
+ sub_domain0_vertices[2],
+ sub_domain0_vertices[3], #
interface12_vertices[0]]
}
# subdomain2
-sub_domain2_vertices = [df.Point(0.0,-1.0),
- df.Point(1.0,-1.0),
+sub_domain2_vertices = [sub_domain0_vertices[0],
+ sub_domain0_vertices[1],
interface12_vertices[1],
interface12_vertices[0] ]
subdomain2_outer_boundary_verts = {
0: [interface12_vertices[0], #
- df.Point(0.0,-1.0), #
- df.Point(1.0,-1.0), #
+ sub_domain0_vertices[0],
+ sub_domain0_vertices[1],
interface12_vertices[1]]
}
+
# subdomain2_outer_boundary_verts = {
# 0: [interface12_vertices[0], df.Point(0.0,0.0)],#
# 1: [df.Point(0.0,0.0), df.Point(1.0,0.0)], #
@@ -85,10 +86,10 @@ isRichards = {
}
##################### MESH #####################################################
-mesh_resolution = 5
+mesh_resolution = 10
##################### TIME #####################################################
timestep_size = 0.01
-number_of_timesteps = 300
+number_of_timesteps = 100
# decide how many timesteps you want analysed. Analysed means, that we write out
# subsequent errors of the L-iteration within the timestep.
number_of_timesteps_to_analyse = 11
@@ -107,6 +108,7 @@ densities = {#
2 : {'wetting' :1}
}
+# gravity_acceleration = 9.81
porosity = {#
# subdom_num : porosity
@@ -122,8 +124,8 @@ L = {#
lambda_param = {#
# subdom_num : lambda parameter for the L-scheme
- 1 : {'wetting' :40},#
- 2 : {'wetting' :40}
+ 1 : {'wetting' :4},#
+ 2 : {'wetting' :4}
}
## relative permeabilty functions on subdomain 1