From 8da2ed26b481c612874ddbc1539d0b177a9bda27 Mon Sep 17 00:00:00 2001
From: David Seus <david.seus@ians.uni-stuttgart.de>
Date: Sat, 28 Sep 2019 18:17:18 +0200
Subject: [PATCH] reimplement gli terms for TPR coupling and run test case
---
LDDsimulation/domainPatch.py | 458 +++++++++---------
.../TP-R-2-patch-test.py | 197 +++++---
2 files changed, 373 insertions(+), 282 deletions(-)
diff --git a/LDDsimulation/domainPatch.py b/LDDsimulation/domainPatch.py
index 6e69ec9..d7bf0ec 100644
--- a/LDDsimulation/domainPatch.py
+++ b/LDDsimulation/domainPatch.py
@@ -344,61 +344,50 @@ class DomainPatch(df.SubDomain):
# 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
+ if self.has_interface is not None:
+ # 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))
+ # if our subdomain assumes Richards and our neighbour TP, we need to
+ # call self.calc_gli_term for the nonwetting phase even though we
+ # are not using it to determin the nonwetting pressure on this patch.
+ # The self.calc_gli_term method saves the calculated values for
+ # communication to the interface dictionries and
+ # the neighbour needs to read this updated gli term,
+ # to approximate the zero pressure.
+ if self.isRichards:
+ # note, that in this case we know that phase == "wetting",
+ # since the solver runs only over phases beloning to our
+ # subdomain. But if self.isRichard, self.has_phases = ["wetting"]
+ neighbour_assumes_R = self.interface[interface].neighbour_isRichards[self.subdomain_index]
+ if not neighbour_assumes_R:
+ unused_nonwetting_gli = self.calc_gli_term(interface_index=interface, phase="nonwetting")
+
+
+ # form2 and rhs2 are different for wetting and nonwetting
# 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":
- if self.has_interface is not None:
- # 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
-
- elif phase == "nonwetting":
- if self.has_interface is not None:
- # 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]:
- # 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)
- 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))
-
+ else:
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, 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
+
if self.has_interface is not None:
form = form1 + form2 + sum(interface_forms)
if self.include_gravity:
@@ -459,71 +448,62 @@ class DomainPatch(df.SubDomain):
# phi_a = self.testfunction[phase]
if self.include_gravity:
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
- # nonwetting phase and assemble the terms as zero form. We don't
- # need to do anything, since gli_term_prev has been
- # initialised to zero.
- pass
- # gli_assembled_tmp[phase] += df.assemble(df.Constant(0)*ds).get_local()
+ if phase == 'wetting':
+ # the wetting phase is always present and we always need
+ # to calculate a gl0 term.
+ # TODO: add intrinsic permeabilty!!!!
+ if self.include_gravity:
+ # z_a included the minus sign, 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:
- if phase == 'wetting':
- # the wetting phase is always present and we always need
- # to calculate a gl0 term.
- # TODO: add intrinsic permeabilty!!!!
- if self.include_gravity:
- # z_a included the minus sign, 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)
+ # phase == 'nonwetting'
+ # we need anothre flux in this case
+ if self.include_gravity:
+ flux = -1/mu_a*ka(1-S(pc_prev))*df.grad(pa_prev - z_a)
else:
- # phase == 'nonwetting'
- # we need anothre flux in this case
- 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)
-
- flux_function = df.project(flux, self.function_space["flux"][phase])
- flux_x, flux_y = flux_function.split()
- gl0 = df.Function(self.function_space["gli"][phase])
- # # get dictionaries of global dof indices and coordinates along
- # # the interface. These dictionaries have the facet indices of
- # # facets belonging to the interface as keys (with respect to
- # # the submesh numbering) and the dictionary
- # # containing pairs {global_dof_index: dof_coordinate} for all
- # # dofs along the facet as values.
- gli_dof_coordinates = self.interface_dof_indices_and_coordinates["gli"][interf_ind][phase]
- # p_dof_coordinates = self.interface_dof_indices_and_coordinates["pressure"][interf_ind][phase]
- # flux_dof_coordinates = self.interface_dof_indices_and_coordinates["flux"][interf_ind][phase]
- # # the attributes local and global for facet numbering mean
- # # the following: local facet index is the index of the facet
- # # with respect to the numbering of the submesh belonging to
- # # the subdomain. global facet index refers to the facet numbering
- # # of the mesh of the whole computational domain.
- local2global_facet = interface.local_to_global_facet_map[subdomain]
- corresponding_dof_index = self.interface_corresponding_dof_index[interf_ind][phase]
- for local_facet_ind, normal in interface.outer_normals[subdomain].items():
- gli_dofs_along_facet = gli_dof_coordinates[local_facet_ind]
- global_facet_ind = local2global_facet[local_facet_ind]
- for gli_dof_index, gli_dof_coord in gli_dofs_along_facet.items():
- flux_x_dof_index = corresponding_dof_index[local_facet_ind][gli_dof_index]["flux_x"]
- flux_y_dof_index = corresponding_dof_index[local_facet_ind][gli_dof_index]["flux_y"]
- p_dof_index = corresponding_dof_index[local_facet_ind][gli_dof_index]["pressure"]
-
- gl0.vector()[gli_dof_index] = flux_x.vector()[flux_x_dof_index]*normal[0] \
- + flux_y.vector()[flux_y_dof_index]*normal[1] \
- - Lambda[phase]*pa_prev.vector()[p_dof_index]
- if debug:
- print(f"gl0.vector()[gli_dof_ind] = {gl0.vector()[gli_dof_index]}")
-
- # save this newly calculated dof to the interface
- # dictionary gli_term_prev
- dof_coord_tupel = (gli_dof_coord[0], gli_dof_coord[1])
- interface.gli_term_prev[subdomain][phase][global_facet_ind].update(
- {dof_coord_tupel: gl0.vector()[gli_dof_index]}
- )
+ flux = -1/mu_a*ka(1-S(pc_prev))*df.grad(pa_prev)
+
+ flux_function = df.project(flux, self.function_space["flux"][phase])
+ flux_x, flux_y = flux_function.split()
+ gl0 = df.Function(self.function_space["gli"][phase])
+ # # get dictionaries of global dof indices and coordinates along
+ # # the interface. These dictionaries have the facet indices of
+ # # facets belonging to the interface as keys (with respect to
+ # # the submesh numbering) and the dictionary
+ # # containing pairs {global_dof_index: dof_coordinate} for all
+ # # dofs along the facet as values.
+ gli_dof_coordinates = self.interface_dof_indices_and_coordinates["gli"][interf_ind][phase]
+ # p_dof_coordinates = self.interface_dof_indices_and_coordinates["pressure"][interf_ind][phase]
+ # flux_dof_coordinates = self.interface_dof_indices_and_coordinates["flux"][interf_ind][phase]
+ # # the attributes local and global for facet numbering mean
+ # # the following: local facet index is the index of the facet
+ # # with respect to the numbering of the submesh belonging to
+ # # the subdomain. global facet index refers to the facet numbering
+ # # of the mesh of the whole computational domain.
+ local2global_facet = interface.local_to_global_facet_map[subdomain]
+ corresponding_dof_index = self.interface_corresponding_dof_index[interf_ind][phase]
+ for local_facet_ind, normal in interface.outer_normals[subdomain].items():
+ gli_dofs_along_facet = gli_dof_coordinates[local_facet_ind]
+ global_facet_ind = local2global_facet[local_facet_ind]
+ for gli_dof_index, gli_dof_coord in gli_dofs_along_facet.items():
+ flux_x_dof_index = corresponding_dof_index[local_facet_ind][gli_dof_index]["flux_x"]
+ flux_y_dof_index = corresponding_dof_index[local_facet_ind][gli_dof_index]["flux_y"]
+ p_dof_index = corresponding_dof_index[local_facet_ind][gli_dof_index]["pressure"]
+
+ gl0.vector()[gli_dof_index] = flux_x.vector()[flux_x_dof_index]*normal[0] \
+ + flux_y.vector()[flux_y_dof_index]*normal[1] \
+ - Lambda[phase]*pa_prev.vector()[p_dof_index]
+ if debug:
+ print(f"gl0.vector()[gli_dof_ind] = {gl0.vector()[gli_dof_index]}")
+
+ # save this newly calculated dof to the interface
+ # dictionary gli_term_prev
+ dof_coord_tupel = (gli_dof_coord[0], gli_dof_coord[1])
+ interface.gli_term_prev[subdomain][phase][global_facet_ind].update(
+ {dof_coord_tupel: gl0.vector()[gli_dof_index]}
+ )
### END calc_gl0_term
@@ -646,115 +626,114 @@ class DomainPatch(df.SubDomain):
# save the iteration number of the neighbour
neighbour_iter_num = interface.current_iteration[neighbour]
- 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 these values 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.")
-
- # dt = self.timestep_size
- Lambda = self.lambda_param[phase]
-
- gli_interface_dofs = self.interface_dof_indices_and_coordinates["gli"][interface_index][phase]
- # p_interface_dofs = self.interface_dof_indices_and_coordinates["pressure"][interface_index][phase]
- local2global_facet = interface.local_to_global_facet_map[subdomain]
- for local_facet_index, gli_facet_dofs2coordinates in gli_interface_dofs.items():
- global_facet_index = local2global_facet[local_facet_index]
- # read gli_prev term from the neighbour.
- gli_prev = interface.gli_term_prev[neighbour][phase][global_facet_index]
- p_prev = pressure_prev[global_facet_index]
- gli_save_dict = gli_save_dictionary[global_facet_index]
- # print(f"\non subdomain{self.subdomain_index}")
- # print(f"on facet with glob_index: {global_facet_index} we have gli_dofs")
- #
- for gli_dof_index, gli_dof_coord in gli_facet_dofs2coordinates.items():
- # find the right previous pressure dof of the interface
- # dictionary to use for the calculation.
- # print(f"coordinates of gli_dof with index {gli_dof_index} I try to match is: {gli_dof_coord}" )
- p_previous_dof = None
- for p_dof_coord_tupel, p_prev_dof in p_prev.items():
- p_coord = np.array([p_dof_coord_tupel[0], p_dof_coord_tupel[1]])
- # print(f"coordinates of p_dof I test this against: {p_coord}")
- a_close_b = np.allclose(gli_dof_coord, p_coord, rtol=1e-10, atol=1e-14)
- b_close_a = np.allclose(p_coord, gli_dof_coord, rtol=1e-10, atol=1e-14)
- # print(f"type(gli_dof_coord) = {type(gli_dof_coord)}")
- # print(f"np.allclose({gli_dof_coord}, {p_coord}, rtol=1e-10, atol=1e-14): {a_close_b }")
- # print(f"np.allclose({p_coord}, {gli_dof_coord}, rtol=1e-10, atol=1e-14): {b_close_a}")
- if a_close_b and b_close_a:
- p_previous_dof = p_prev_dof
-
- if p_previous_dof is None:
- raise RuntimeError(f"didn't find p_previous_dof corresponding to dict key ({gli_dof_coord})",
- "something is wrong")
- # same with the gli_previous_dof
- gli_previous_dof = None
- for gli_prev_coord_tupel, gli_prev_dof in gli_prev.items():
- gli_prev_coord = np.array([gli_prev_coord_tupel[0], gli_prev_coord_tupel[1]])
- # due to the warning in the documentation of np.allclose
- # that np.allclose is not symetric, we test if both
- # np.allclose(a,b) and np.allclose(b,a) are true.
- a_close_b = np.allclose(gli_dof_coord, gli_prev_coord, rtol=1e-10, atol=1e-14)
- b_close_a = np.allclose(gli_prev_coord, gli_dof_coord, rtol=1e-10, atol=1e-14)
- if a_close_b and b_close_a:
- gli_previous_dof = gli_prev_dof
- if gli_previous_dof is None:
- raise RuntimeError(f"didn't find gli_previous_dof corresponding to dict key ({gli_dof_coord})",
- "something is wrong")
-
- gli_value = -2*Lambda*p_previous_dof - gli_previous_dof
- gli.vector()[gli_dof_index] = gli_value
- gli_dof_coord_tupel = (gli_dof_coord[0], gli_dof_coord[1])
- gli_save_dict.update({gli_dof_coord_tupel: 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][phase].update(#
- {global_facet_index: interface.gli_term[neighbour][phase][global_facet_index].copy()}
- )
- interface.pressure_values_prev[neighbour][phase].update(
- {global_facet_index: interface.pressure_values[neighbour][phase][global_facet_index].copy()}
+ # 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 these values 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.")
+
+ # dt = self.timestep_size
+ Lambda = self.lambda_param[phase]
+
+ gli_interface_dofs = self.interface_dof_indices_and_coordinates["gli"][interface_index][phase]
+ # p_interface_dofs = self.interface_dof_indices_and_coordinates["pressure"][interface_index][phase]
+ local2global_facet = interface.local_to_global_facet_map[subdomain]
+ for local_facet_index, gli_facet_dofs2coordinates in gli_interface_dofs.items():
+ global_facet_index = local2global_facet[local_facet_index]
+ # read gli_prev term from the neighbour.
+ gli_prev = interface.gli_term_prev[neighbour][phase][global_facet_index]
+ p_prev = pressure_prev[global_facet_index]
+ gli_save_dict = gli_save_dictionary[global_facet_index]
+ # print(f"\non subdomain{self.subdomain_index}")
+ # print(f"on facet with glob_index: {global_facet_index} we have gli_dofs")
+ #
+ for gli_dof_index, gli_dof_coord in gli_facet_dofs2coordinates.items():
+ # find the right previous pressure dof of the interface
+ # dictionary to use for the calculation.
+ # print(f"coordinates of gli_dof with index {gli_dof_index} I try to match is: {gli_dof_coord}" )
+ p_previous_dof = None
+ for p_dof_coord_tupel, p_prev_dof in p_prev.items():
+ p_coord = np.array([p_dof_coord_tupel[0], p_dof_coord_tupel[1]])
+ # print(f"coordinates of p_dof I test this against: {p_coord}")
+ a_close_b = np.allclose(gli_dof_coord, p_coord, rtol=1e-10, atol=1e-14)
+ b_close_a = np.allclose(p_coord, gli_dof_coord, rtol=1e-10, atol=1e-14)
+ # print(f"type(gli_dof_coord) = {type(gli_dof_coord)}")
+ # print(f"np.allclose({gli_dof_coord}, {p_coord}, rtol=1e-10, atol=1e-14): {a_close_b }")
+ # print(f"np.allclose({p_coord}, {gli_dof_coord}, rtol=1e-10, atol=1e-14): {b_close_a}")
+ if a_close_b and b_close_a:
+ p_previous_dof = p_prev_dof
+
+ if p_previous_dof is None:
+ raise RuntimeError(f"didn't find p_previous_dof corresponding to dict key ({gli_dof_coord})",
+ "something is wrong")
+ # same with the gli_previous_dof
+ gli_previous_dof = None
+ for gli_prev_coord_tupel, gli_prev_dof in gli_prev.items():
+ gli_prev_coord = np.array([gli_prev_coord_tupel[0], gli_prev_coord_tupel[1]])
+ # due to the warning in the documentation of np.allclose
+ # that np.allclose is not symetric, we test if both
+ # np.allclose(a,b) and np.allclose(b,a) are true.
+ a_close_b = np.allclose(gli_dof_coord, gli_prev_coord, rtol=1e-10, atol=1e-14)
+ b_close_a = np.allclose(gli_prev_coord, gli_dof_coord, rtol=1e-10, atol=1e-14)
+ if a_close_b and b_close_a:
+ gli_previous_dof = gli_prev_dof
+ if gli_previous_dof is None:
+ raise RuntimeError(f"didn't find gli_previous_dof corresponding to dict key ({gli_dof_coord})",
+ "something is wrong")
+
+ gli_value = -2*Lambda*p_previous_dof - gli_previous_dof
+ gli.vector()[gli_dof_index] = gli_value
+ gli_dof_coord_tupel = (gli_dof_coord[0], gli_dof_coord[1])
+ gli_save_dict.update({gli_dof_coord_tupel: 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][phase].update(#
+ {global_facet_index: interface.gli_term[neighbour][phase][global_facet_index].copy()}
)
+ interface.pressure_values_prev[neighbour][phase].update(
+ {global_facet_index: interface.pressure_values[neighbour][phase][global_facet_index].copy()}
+ )
return gli
### END calc_gli_term
@@ -799,7 +778,23 @@ class DomainPatch(df.SubDomain):
pressure_space = df.FunctionSpace(self.mesh, 'P', degree)
gli_space = df.FunctionSpace(self.mesh, 'DG', degree)
flux_space = df.VectorFunctionSpace(self.mesh, 'DG', degree)
- for phase in self.has_phases:
+
+ if self.isRichards:
+ # if our domain assumes Richards but the neighbour of the interface
+ # assumes TP, we need to initialise interface values for the nonwetting
+ # phase, too, for communication, see the article. Therefor we need to
+ # initialise the nonwetting function space too.
+ for interf in self.has_interface:
+ neighbour_assumes_R = self.interface[interf].neighbour_isRichards[self.subdomain_index]
+ if not neighbour_assumes_R:
+ subdom_has_phases = ["wetting", "nonwetting"]
+ # if only one neighbour assumes TP we ste up the nonwetting
+ # space.
+ break
+ else:
+ subdom_has_phases = self.has_phases
+
+ for phase in subdom_has_phases:
if function == "pressure":
self.function_space[function].update({phase: pressure_space})
self.trialfunction.update({phase: df.TrialFunction(self.function_space["pressure"][phase])})
@@ -924,7 +919,16 @@ class DomainPatch(df.SubDomain):
)
if debug:
print(f"\nOn subdomain {self.subdomain_index}")
- for phase in self.has_phases:
+
+ # if our subdomain is Richards but the neighbour is TP, we
+ # need to calculate interface dof indices for communication
+ # also for the nonwetting phase
+ if self.isRichards and not interface.neighbour_isRichards[self.subdomain_index]:
+ subdom_has_phases = ["wetting", "nonwetting"]
+ else:
+ subdom_has_phases = self.has_phases
+
+ for phase in subdom_has_phases:
if function == "flux":
self.interface_dof_indices_and_coordinates[function][interface_index].update(#
{phase: dict()}
@@ -957,7 +961,16 @@ class DomainPatch(df.SubDomain):
self.interface_corresponding_dof_index.update(
{interf_ind: dict()}
)
- for phase in self.has_phases:
+
+ # if our subdomain is Richards but the neighbour is TP, we
+ # need to calculate interface dof indices for communication
+ # also for the nonwetting phase
+ if self.isRichards and not interface.neighbour_isRichards[self.subdomain_index]:
+ subdom_has_phases = ["wetting", "nonwetting"]
+ else:
+ subdom_has_phases = self.has_phases
+
+ for phase in subdom_has_phases:
self.interface_corresponding_dof_index[interf_ind].update(
{phase: dict()}
)
@@ -1098,15 +1111,24 @@ class DomainPatch(df.SubDomain):
for interf in self.has_interface:
marker_value = self.interface[interf].marker_value
space = self.function_space
+
+ neighbour_assumes_R = self.interface[interf].neighbour_isRichards[self.subdomain_index]
+ # if our domain assumes Richards but the neighbour of the interface
+ # assumes TP, we need to initialise interface values for the nonwetting
+ # phase, too for communication, see the article.
+ if self.isRichards and not neighbour_assumes_R:
+ subdom_has_phases=["wetting","nonwetting"]
+ else:
+ subdom_has_phases=self.has_phases
+
self.interface[interf].init_interface_values(
interface_marker=marker,
interface_marker_value=marker_value,
function_space=space,
subdomain_index=self.subdomain_index,
- has_phases=self.has_phases,
+ has_phases=subdom_has_phases,
)
-
def _calc_gravity_expressions(self):
""" calculate the gravity term if it is included"""
g = df.Constant(self.gravity_acceleration) #
diff --git a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-test.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-test.py
index 62645f7..cd8315d 100755
--- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-test.py
+++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-test.py
@@ -8,34 +8,78 @@ import domainPatch as dp
import LDDsimulation as ldd
import functools as ft
import helpers as hlp
+import datetime
+import os
+import pandas as pd
+
+date = datetime.datetime.now()
+datestr = date.strftime("%Y-%m-%d")
#import ufl as ufl
# init sympy session
sym.init_printing()
-use_case = "TP-R-two-patch-realistic "
-solver_tol = 1E-6
-
-############ GRID #######################ΓΌ
-mesh_resolution = 20
-timestep_size = 0.00001
-number_of_timesteps = 15
+use_case = "TP-R-2-patch-realistic-new-implementation"
+# solver_tol = 6E-7
+max_iter_num = 1000
+FEM_Lagrange_degree = 1
+mesh_study = False
+resolutions = {
+ # 1: 7e-7,
+ # 2: 7e-7,
+ # 4: 7e-7,
+ # 8: 7e-7,
+ # 16: 7e-7,
+ 32: 7e-7,
+ # 64: 7e-7,
+ # 128: 7e-7,
+ # 256: 7e-7
+ }
+
+############ GRID #######################
+# mesh_resolution = 20
+timestep_size = 0.001
+number_of_timesteps = 20
+plot_timestep_every = 1
# 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 = 5
-starttime = 0
+number_of_timesteps_to_analyse = 10
+starttimes = [0.0]
-Lw = 1 #/timestep_size
+Lw = 0.025 #/timestep_size
Lnw=Lw
-l_param_w = 40
-l_param_nw = 40
+lambda_w = 40
+lambda_nw = 40
-include_gravity = True
+include_gravity = False
debugflag = True
-analyse_condition = False
+analyse_condition = True
+
+output_string = "./output/{}-{}_timesteps{}_P{}".format(datestr, use_case, number_of_timesteps, FEM_Lagrange_degree)
+
+# toggle what should be written to files
+if mesh_study:
+ write_to_file = {
+ 'space_errornorms': True,
+ 'meshes_and_markers': True,
+ 'L_iterations_per_timestep': False,
+ 'solutions': False,
+ 'absolute_differences': False,
+ 'condition_numbers': analyse_condition,
+ 'subsequent_errors': False
+ }
+else:
+ write_to_file = {
+ 'space_errornorms': True,
+ 'meshes_and_markers': True,
+ 'L_iterations_per_timestep': False,
+ 'solutions': True,
+ 'absolute_differences': True,
+ 'condition_numbers': analyse_condition,
+ 'subsequent_errors': True
+ }
-output_string = "./output/after_reimplementing_gravity_term-realistic_number_of_timesteps{}_".format(number_of_timesteps)
##### Domain and Interface ####
# global simulation domain domain
@@ -139,10 +183,10 @@ L = {#
lambda_param = {#
# subdom_num : lambda parameter for the L-scheme
- 1 : {'wetting' :l_param_w},
- # 'nonwetting': l_param},#
- 2 : {'wetting' :l_param_w,
- 'nonwetting': l_param_nw}
+ 1 : {'wetting' :lambda_w,
+ 'nonwetting': lambda_nw},#
+ 2 : {'wetting' :lambda_w,
+ 'nonwetting': lambda_nw}
}
## relative permeabilty functions on subdomain 1
@@ -325,7 +369,7 @@ t = sym.symbols('t', positive=True)
p_e_sym = {
1: {'wetting': (-5.0 - (1.0 + t*t)*(1.0 + x*x + y*y))}, #*(1-x)**2*(1+x)**2*(1-y)**2},
2: {'wetting': (-5.0 - (1.0 + t*t)*(1.0 + x*x)), #*(1-x)**2*(1+x)**2*(1+y)**2,
- 'nonwetting': (-1-t*(1.1+y + x**2))*y**2}, #*(1-x)**2*(1+x)**2*(1+y)**2},
+ 'nonwetting': (-1-t*(1.1+y + x**2))*y**3}, #*(1-x)**2*(1+x)**2*(1+y)**2},
} #-y*y*(sym.sin(-2*t+2*x)*sym.sin(1/2*y-1.2*t)) - t*t*x*(0.5-y)*y*(1-x)
@@ -391,47 +435,72 @@ for subdomain in isRichards.keys():
)
-# def saturation(pressure, subdomain_index):
-# # inverse capillary pressure-saturation-relationship
-# return df.conditional(pressure < 0, 1/((1 - pressure)**(1/(subdomain_index + 1))), 1)
-#
-# sa
-
-write_to_file = {
- 'meshes_and_markers': True,
- 'L_iterations': True
-}
-
-
-# initialise LDD simulation class
-simulation = ldd.LDDsimulation(tol = 1E-14, LDDsolver_tol=solver_tol, debug=debugflag)
-simulation.set_parameters(use_case=use_case,
- output_dir = output_string,#
- subdomain_def_points = subdomain_def_points,#
- isRichards = isRichards,#
- interface_def_points = interface_def_points,#
- outer_boundary_def_points = outer_boundary_def_points,#
- adjacent_subdomains = adjacent_subdomains,#
- mesh_resolution = mesh_resolution,#
- viscosity = viscosity,#
- porosity = porosity,#
- L = L,#
- lambda_param = lambda_param,#
- relative_permeability = relative_permeability,#
- saturation = sat_pressure_relationship,#
- starttime = starttime,#
- number_of_timesteps = number_of_timesteps,
- number_of_timesteps_to_analyse = number_of_timesteps_to_analyse,
- timestep_size = timestep_size,#
- sources = source_expression,#
- initial_conditions = initial_condition,#
- dirichletBC_expression_strings = dirichletBC,#
- exact_solution = exact_solution,#
- densities=densities,
- include_gravity=include_gravity,
- write2file = write_to_file,#
- )
-
-simulation.initialise()
-# simulation.write_exact_solution_to_xdmf()
-simulation.run(analyse_condition=analyse_condition)
+for starttime in starttimes:
+ for mesh_resolution, solver_tol in resolutions.items():
+ # initialise LDD simulation class
+ simulation = ldd.LDDsimulation(
+ tol=1E-14,
+ LDDsolver_tol=solver_tol,
+ debug=debugflag,
+ max_iter_num=max_iter_num,
+ FEM_Lagrange_degree=FEM_Lagrange_degree,
+ mesh_study=mesh_study
+ )
+
+ simulation.set_parameters(use_case=use_case,
+ output_dir=output_string,
+ subdomain_def_points=subdomain_def_points,
+ isRichards=isRichards,
+ interface_def_points=interface_def_points,
+ outer_boundary_def_points=outer_boundary_def_points,
+ adjacent_subdomains=adjacent_subdomains,
+ mesh_resolution=mesh_resolution,
+ viscosity=viscosity,
+ porosity=porosity,
+ L=L,
+ lambda_param=lambda_param,
+ relative_permeability=relative_permeability,
+ saturation=sat_pressure_relationship,
+ starttime=starttime,
+ number_of_timesteps=number_of_timesteps,
+ number_of_timesteps_to_analyse=number_of_timesteps_to_analyse,
+ plot_timestep_every=plot_timestep_every,
+ timestep_size=timestep_size,
+ sources=source_expression,
+ initial_conditions=initial_condition,
+ dirichletBC_expression_strings=dirichletBC,
+ exact_solution=exact_solution,
+ densities=densities,
+ include_gravity=include_gravity,
+ write2file=write_to_file,
+ )
+
+ simulation.initialise()
+ output_dir = simulation.output_dir
+ # simulation.write_exact_solution_to_xdmf()
+ output = simulation.run(analyse_condition=analyse_condition)
+ for subdomain_index, subdomain_output in output.items():
+ mesh_h = subdomain_output['mesh_size']
+ for phase, different_errornorms in subdomain_output['errornorm'].items():
+ filename = output_dir + "subdomain{}-space-time-errornorm-{}-phase.csv".format(subdomain_index, phase)
+ # for errortype, errornorm in different_errornorms.items():
+
+ # eocfile = open("eoc_filename", "a")
+ # eocfile.write( str(mesh_h) + " " + str(errornorm) + "\n" )
+ # eocfile.close()
+ # if subdomain.isRichards:mesh_h
+ data_dict = {
+ 'mesh_parameter': mesh_resolution,
+ 'mesh_h': mesh_h,
+ }
+ for error_type, errornorms in different_errornorms.items():
+ data_dict.update(
+ {error_type: errornorms}
+ )
+ errors = pd.DataFrame(data_dict, index=[mesh_resolution])
+ # check if file exists
+ if os.path.isfile(filename) == True:
+ with open(filename, 'a') as f:
+ errors.to_csv(f, header=False, sep='\t', encoding='utf-8', index=False)
+ else:
+ errors.to_csv(filename, sep='\t', encoding='utf-8', index=False)
--
GitLab