diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index 68aa1a00ecd70ca23c25265b6fd8418dc052b469..a076bc744e8b1a1000a560fa75a25d1feaa07569 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -304,7 +304,7 @@ class LDDsimulation(object):
self.timestep_num = int(0)
# note that at this point of the code self.t = self.starttime as set in
# the beginning.
- while self.timestep_num < self.number_of_timesteps:
+ while self.timestep_num <= self.number_of_timesteps:
print(f"entering timestep ",
"{}".format(self.timestep_num),
"at time t = "+ "{number:.{digits}f}".format(number = self.t, digits = 4))
@@ -671,8 +671,8 @@ class LDDsimulation(object):
evaluate the exact solution of the simulation (in case there is one) at
all timesteps and write it to the solution file.
"""
- t = copy.copy(self.starttime)
for subdom_ind, subdomain in self.subdomain.items():
+ t = copy.copy(self.starttime)
file = self.solution_file[subdom_ind]
for timestep in range(0,self.number_of_timesteps+1):
for phase in subdomain.has_phases:
@@ -683,7 +683,7 @@ class LDDsimulation(object):
file.write(tmp_exact_pressure, t)
t += self.timestep_size
- file.close()
+
def write_subsequent_errors_to_csv(self,#
filename: str, #
diff --git a/LDDsimulation/boundary_and_interface.py b/LDDsimulation/boundary_and_interface.py
index 57e5a220d4d962f17ace7e41dad9dac042bf6209..b50d921cf6894381ed01377a3ba78f78de225f71 100644
--- a/LDDsimulation/boundary_and_interface.py
+++ b/LDDsimulation/boundary_and_interface.py
@@ -562,216 +562,6 @@ class Interface(BoundaryPart):
# end init_interface_values
- # def read_pressure_dofs(self, from_function: df.Function, #
- # interface_dofs: np.array,#
- # dof_to_vert_map: np.array,#
- # local_to_parent_vertex_map: np.array,#
- # phase: str,#
- # subdomain_ind: int,#
- # previous_iter: bool = False,#
- # ) -> None:
- # """ read dofs of from_function with indices interface_dofs and write to self.pressure_values
- #
- # Read the degrees of freedom of the function from_function with indices
- # interface_dofs, i.e. the dofs of the interface, and write them to one of the dictionaries
- # of the interface
- # self.pressure_values[subdomain_ind][phase], (if pressure == True)
- # self.pressure_values_prev[subdomain_ind][phase], (if previous_iter == True)
- #
- # # Parameters
- # from_function: #type: df.Function, function to save the dofs of
- # interface_dofs: #type: np.array, index array of dofs to be saved
- # dof_to_vert_map: #type: np.array, dof to vert map of the function
- # space on the submesh.
- # local_to_parent_vertex_map: #type: np.array,
- # phase #type: str is either 'wetting' or
- # 'nonwetting' and determins
- # the dict to be written to.
- # subdomain_ind #type: int subdomain index
- # previous_iter #type: bool determins wether a previous
- # iterate is being written or not
- # """
- # if self.pressure_values == None:
- # raise RuntimeError("self.pressure_values not initiated. You forgot to \
- # run self.init_interface_values")
- #
- # parent = local_to_parent_vertex_map
- # d2v = dof_to_vert_map
- # if not previous_iter:
- # for dof_index in interface_dofs:
- # # from_function.vector() 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.
- # self.pressure_values[subdomain_ind][phase].update({parent[d2v[dof_index]]: from_function.vector()[dof_index]})
- # # print(f"have written subdomain{subdomain_ind} pressure to interface{self.global_index} for phase {phase}\n", self.pressure_values[subdomain_ind][phase])
- # else:
- # for dof_index in interface_dofs:
- # # from_function.vector() 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.
- # self.pressure_values_prev[subdomain_ind][phase].update({parent[d2v[dof_index]]: from_function.vector()[dof_index]})
- # # print(f"have written subdomain{subdomain_ind} pressure_prev_iter to interface{self.global_index} for phase {phase}\n", self.pressure_values_prev[subdomain_ind][phase])
- #
- #
- # def read_gli_dofs(self, from_array: np.array, #
- # interface_dofs: np.array,#
- # dof_to_vert_map: np.array,#
- # local_to_parent_vertex_map: np.array,#
- # phase: str,#
- # subdomain_ind: int,#
- # previous_iter: bool = False,#
- # ) -> None:
- # """ Read dofs of from_array with indices interface_dofs and write to self.gli_term
- #
- # Read the degrees of freedom of the array from_array with indices
- # interface_dofs, i.e. the dofs of the interface, and write them to one of the dictionaries
- # of the interface
- # self.gli_term[subdomain_ind][phase], (if gl == True)
- # self.gli_term_prev[subdomain_ind][phase] (if previous_iter == True)
- #
- # # Parameters
- # from_array: #type: np.array, array to save the dofs of
- # interface_dofs: #type: np.array, index array of dofs to be saved
- # dof_to_vert_map: #type: np.array, dof to vert map of the function
- # space on the submesh.
- # local_to_parent_vertex_map: #type: np.array,
- # phase #type: str is either 'wetting' or
- # 'nonwetting' and determins
- # the dict to be written to.
- # subdomain_ind #type: int subdomain index
- # previous_iter #type: bool determins wether a previous
- # iterate is being written or not
- # """
- # if self.pressure_values == None:
- # raise RuntimeError("self.pressure_values not initiated. You forgot to \
- # run self.init_interface_values")
- #
- # parent = local_to_parent_vertex_map
- # d2v = dof_to_vert_map
- # if not previous_iter:
- # for dof_index in interface_dofs:
- # # 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.
- # self.gli_term[subdomain_ind][phase].update({parent[d2v[dof_index]]: from_array[dof_index]})
- # # print(f"have written subdomain{subdomain_ind} gl dofs to interface{self.global_index} for phase {phase}\n", self.gli_term[subdomain_ind][phase])
- # else:
- # for dof_index in interface_dofs:
- # # 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.
- # self.gli_term_prev[subdomain_ind][phase].update({parent[d2v[dof_index]]: from_array[dof_index]})
- # # print(f"have written subdomain{subdomain_ind} gl_prev dofs to interface{self.global_index} for phase {phase}\n", self.gli_term_prev[subdomain_ind][phase])
- #
- #
- # def write_pressure_dofs(self, to_function: df.Function, #
- # interface_dofs: np.array,#
- # dof_to_vert_map: np.array,#
- # local_to_parent_vertex_map: np.array,#
- # phase: str,#
- # subdomain_ind: int,#
- # # pressure: bool = False,#
- # # gl: bool = False,#
- # previous_iter: bool = False,#
- # ) -> None:
- # """ read dofs off self.pressure_values and overwrite the dofs of to_function
- # with indices interface_dofs
- #
- # Read the degrees of freedom of one of the dictionaries of the interface
- # self.pressure_values[subdomain_ind][phase], (if pressure == True)
- # self.pressure_values_prev[subdomain_ind][phase], (if previous_iter == True)
- # and overwrite the dofs of the function to_function with indices
- # interface_dofs, i.e. update the interface dofs of the function.
- # .
- #
- # # Parameters
- # to_function: #type: df.Function, function to overwrite the
- # dofs of.
- # interface_dofs: #type: np.array, index array of interface
- # dofs of from_function
- # dof_to_vert_map: #type: np.array, dof to vert map of the function
- # space on the submesh.
- # local_to_parent_vertex_map: #type: np.array,
- # phase #type: str is either 'wetting' or
- # 'nonwetting' and determins
- # the dict to be read of.
- # subdomain_ind #type: int subdomain index
- # previous_iter #type: bool determins wether a previous
- # iterate is being read or not
- # """
- # if self.pressure_values == None:
- # raise RuntimeError("self.pressure_values not initiated. You forgot to \
- # run self.init_interface_values")
- #
- # parent = local_to_parent_vertex_map
- # d2v = dof_to_vert_map
- #
- # if not previous_iter:
- # for dof_index in interface_dofs:
- # # from_function.vector() 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.
- # to_function.vector()[dof_index] = self.pressure_values[subdomain_ind][phase][parent[d2v[dof_index]]]
- # # print("read pressure interface values\n", self.pressure_values[subdomain_ind][phase])
- # # print("wrote these dofs to function \n", to_function.vector()[:])
- # else:
- # for dof_index in interface_dofs:
- # # from_function.vector() 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.
- # to_function.vector()[dof_index] = self.pressure_values_prev[subdomain_ind][phase][parent[d2v[dof_index]]]
- # # print("read previous pressure interface values\n", self.pressure_values_prev[subdomain_ind][phase])
- # # print("wrote these dofs to function \n", to_function.vector()[:])
- #
- #
- # def write_gli_dofs(self, to_array: np.array, #
- # interface_dofs: np.array,#
- # dof_to_vert_map: np.array,#
- # local_to_parent_vertex_map: np.array,#
- # phase: str,#
- # subdomain_ind: int,#
- # previous_iter: bool = False,#
- # ) -> None:
- # """ read dofs off self.gli_term and overwrite the dofs stored in the array
- # to_array with indices interface_dofs
- #
- # Read the degrees of freedom of one of the dictionaries of the interface
- # self.gli_term[subdomain_ind][phase], (if gl == True)
- # self.gli_term_prev[subdomain_ind][phase] (if previous_iter == True)
- # and overwrite the values (holding dofs) of the array to_array with indices
- # interface_dofs, i.e. update the interface dofs of the gli terms.
- # .
- #
- # # Parameters
- # to_array: #type: np.array, array to overwrite the
- # dofs of.
- # interface_dofs: #type: np.array, index array of interface
- # dofs of from_function
- # dof_to_vert_map: #type: np.array, dof to vert map of the function
- # space on the submesh.
- # local_to_parent_vertex_map: #type: np.array,
- # phase #type: str is either 'wetting' or
- # 'nonwetting' and determins
- # the dict to be read of.
- # subdomain_ind #type: int subdomain index
- #
- # previous_iter #type: bool determins wether a previous
- # iterate is being read or not
- # """
- # if self.pressure_values == None:
- # raise RuntimeError("self.pressure_values not initiated. You forgot to \
- # run self.init_interface_values")
- #
- # parent = local_to_parent_vertex_map
- # d2v = dof_to_vert_map
- # if not previous_iter:
- # for dof_index in interface_dofs:
- # # to_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.
- # to_array[dof_index] = self.gli_term[subdomain_ind][phase][parent[d2v[dof_index]]]
- # # print("read gl interface values\n", self.gli_term[subdomain_ind][phase])
- # # print("wrote these dofs to array \n", to_array[:])
- # else:
- # for dof_index in interface_dofs:
- # # to_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.
- # to_array[dof_index] = self.gli_term_prev[subdomain_ind][phase][parent[d2v[dof_index]]]
- # # print("read gl interface values\n", self.gli_term_prev[subdomain_ind][phase])
- # # print("wrote these dofs to array \n", to_array[:])
#### Private Methods ####
diff --git a/LDDsimulation/domainPatch.py b/LDDsimulation/domainPatch.py
index bff23f732bfa5cee2b10b135de22d94cbe1fbe9f..acdb82f3b32807857c967bc2adb2ec9a411b16d6 100644
--- a/LDDsimulation/domainPatch.py
+++ b/LDDsimulation/domainPatch.py
@@ -278,7 +278,6 @@ class DomainPatch(df.SubDomain):
for phase in self.has_phases:
p_dof_coordinates = self.interface_dof_indices_and_coordinates["pressure"][ind][phase]
local2global_facet = interface.local_to_global_facet_map[subdomain]
- # for facet_ind in interface.facets[subdomain].keys():
for subdomain_facet_index, facet_dof_coord_dict in p_dof_coordinates.items():
global_facet_index = local2global_facet[subdomain_facet_index]
# set dictionary to save the values to.
@@ -446,7 +445,7 @@ 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. We don't
- # need to do anything, since gli_assembled_tmp[phase] has been
+ # 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()
@@ -820,7 +819,6 @@ class DomainPatch(df.SubDomain):
interface values for communication.
"""
# we don't need to communicate flux values directly.
- function_name = ["pressure", "gli"]
marker = self.interface_marker
for interf in self.has_interface:
@@ -835,104 +833,6 @@ class DomainPatch(df.SubDomain):
)
- # def _calc_dof_indices_of_interfaces(self):
- # """ calculate dof indices of each interface """
- # V = self.function_space["pressure"]
- # marker = self.interface_marker
- # for ind in self.has_interface:
- # # the marker value for the interfaces is always global interface index+1
- # marker_value = self.interface[ind].marker_value
- # self._dof_indices_of_interface.update(#
- # {ind: dict()}
- # )
- # for phase in self.has_phases:
- # self._dof_indices_of_interface[ind].update(#
- # {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])
- #
- # def _calc_interface_has_common_dof_indices(self):
- # """ populate dictionary self._interface_has_common_dof_indices
- #
- # Determin for each interface in self.has_interface the dof indices that
- # also belong to another interface and store them in the dictionary
- # self._interface_has_common_dof_indices
- # This method populates the dictionary
- # self._interface_has_common_dof_indices
- # """
- # for interf_ind, interf_dofs in self._dof_indices_of_interface.items():
- # # initialise bool numpy arry for the search.
- # if self.isRichards:
- # is_part_of_neighbour_interface_w = np.zeros(interf_dofs['wetting'].size, dtype=bool)
- # else:
- # is_part_of_neighbour_interface_w = np.zeros(interf_dofs['wetting'].size, dtype=bool)
- # is_part_of_neighbour_interface_nw = np.zeros(interf_dofs['nonwetting'].size, dtype=bool)
- #
- # for other_interf_ind, other_interf_dofs in self._dof_indices_of_interface.items():
- # if other_interf_ind is not interf_ind:
- # if self.isRichards:
- # is_part_of_neighbour_interface_w += np.isin(interf_dofs['wetting'], other_interf_dofs['wetting'], assume_unique=True)
- # else:
- # is_part_of_neighbour_interface_w += np.isin(interf_dofs['wetting'], other_interf_dofs['wetting'], assume_unique=True)
- # is_part_of_neighbour_interface_nw += np.isin(interf_dofs['nonwetting'], other_interf_dofs['nonwetting'], assume_unique=True)
- # # this carries a numpyarray of the indices which belong to other
- # # interfaces, too
- # self._interface_has_common_dof_indices.update({interf_ind: dict()})
- # if self.isRichards:
- # self._interface_has_common_dof_indices[interf_ind].update(
- # {'wetting': interf_dofs['wetting'][is_part_of_neighbour_interface_w]}
- # )
- # else:
- # self._interface_has_common_dof_indices[interf_ind].update(
- # {'wetting': interf_dofs['wetting'][is_part_of_neighbour_interface_w],
- # 'nonwetting': interf_dofs['nonwetting'][is_part_of_neighbour_interface_nw]}
- # )
- # # for phase in self.has_phases:
- # # print(f"self._interface_has_common_dof_indices[{interf_ind}][{phase}]", self._interface_has_common_dof_indices[interf_ind][phase])
- # # print(f"python id of self._interface_has_common_dof_indices[{interf_ind}][{phase}]", id(self._interface_has_common_dof_indices[interf_ind][phase]))
-
-
- # def _calc_normal_vectors_norms_for_common_dofs(self):
- # """ calculate norm(n1+n2) where n1 and n2 are normal vectors of adjacent
- # interfaces. This is used by the functions calc_gli_term and calc_gl0_term
- # """
- # self._normal_vector_norms_for_common_dofs = dict()
- # self.mesh.init(0)
- # for interface_ind in self.has_interface:
- # # if one phase has common dofs to two interfaces, the both phases
- # # have, so it suffices to check for the wetting phase as it is always
- # # present.
- # if self._interface_has_common_dof_indices[interface_ind]['wetting'].size != 0:
- # self._normal_vector_norms_for_common_dofs.update({interface_ind: dict()})
- # for phase in self.has_phases:
- # self._normal_vector_norms_for_common_dofs[interface_ind].update(
- # {phase: dict()}
- # )
- # for dof in self._interface_has_common_dof_indices[interface_ind][phase]:
- # common_vertex_index = self.dof2vertex[phase][dof]
- # for index, vertex in enumerate(df.vertices(self.mesh)):
- # if index == common_vertex_index:
- # common_vertex = vertex
- # # print(f"\non subdomain{self.subdomain_index}, interface{interface_ind} and phase {phase}")
- # # print(f"vertex_{index} = {vertex.point()[:]}, type(vertex) = {type(vertex)}")
- # linear_combination_of_normals = np.zeros(np.size(vertex.point()[:]), dtype=float)
- # # normals = []
- # # print(f"common_vertex = {common_vertex}",
- # # f"common_vertex_index = {common_vertex_index}")
- # for facet_containing_vertex in df.facets(common_vertex):#df.cells(common_vertex):
- # if facet_containing_vertex.exterior():
- # # normals.append(facet_containing_vertex.normal()[:])
- # # print(f"type(facet_normal) = {type(facet_containing_vertex.normal()[:])}")
- # linear_combination_of_normals += facet_containing_vertex.normal()[:]
- #
- # # print(f"normals of common_vertex: {common_vertex} are ", normals)
- # # print(f"linear combination is: {linear_combination_of_normals[:]}.")
- # # print(f"with length {np.linalg.norm(linear_combination_of_normals, ord=2, axis=None)}. For comparison, np.sqrt(2) = {np.sqrt(2)}.")
- # self._normal_vector_norms_for_common_dofs[interface_ind][phase].update(
- # {dof: np.linalg.norm(linear_combination_of_normals)}
- # )
-
-
def _calc_gravity_expressions(self):
""" calculate the gravity term if it is included"""
g = df.Constant(self.gravity_acceleration) #
@@ -947,7 +847,6 @@ class DomainPatch(df.SubDomain):
def write_solution_to_xdmf(self, file: tp.Type[SolutionFile], #
time: float = None,#
write_iter_for_fixed_time: bool = False,#
- exact_solution: bool = False,#
):
"""
Weither write the solution of the simulation and corresponding time
@@ -971,11 +870,6 @@ class DomainPatch(df.SubDomain):
if not write_iter_for_fixed_time:
self.pressure[phase].rename("pressure_"+"{}".format(phase), "pressure_"+"{}".format(phase))
file.write(self.pressure[phase], t)
- # if self.pressure_exact:
- # self.pressure_exact[phase].t = t
- # tmp_exact_pressure = df.interpolate(self.pressure_exact[phase], self.function_space["pressure"][phase])
- # tmp_exact_pressure.rename("exact_pressure_"+"{}".format(phase), "exactpressure_"+"{}".format(phase))
- # file.write(tmp_exact_pressure, t)
else:
i = self.iteration_number
self.pressure[phase].rename("pressure_"+"{}".format(phase), #
diff --git a/RR-multi-patch-plus-gravity/RR-multi-patch-with-gravity.py b/RR-multi-patch-plus-gravity/RR-multi-patch-with-gravity.py
index 78d7cfcc272d1c3eb5b0f6af180cf101d9de67d4..b64f006c76508861ba40dc5a882520dc0cdf9324 100755
--- a/RR-multi-patch-plus-gravity/RR-multi-patch-with-gravity.py
+++ b/RR-multi-patch-plus-gravity/RR-multi-patch-with-gravity.py
@@ -15,12 +15,12 @@ sym.init_printing()
# ----------------------------------------------------------------------------#
# ------------------- MESH ---------------------------------------------------#
# ----------------------------------------------------------------------------#
-mesh_resolution = 50
+mesh_resolution = 51
# ----------------------------------------:-------------------------------------#
# ------------------- TIME ---------------------------------------------------#
# ----------------------------------------------------------------------------#
timestep_size = 0.01
-number_of_timesteps = 1000
+number_of_timesteps = 160
# 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
@@ -169,7 +169,7 @@ L = {
4: {'wetting': 0.25}
}
-lamdal_w = 30
+lamdal_w = 40
# subdom_num : lambda parameter for the L-scheme
lambda_param = {
1: {'wetting': lamdal_w},
@@ -435,7 +435,7 @@ write_to_file = {
}
# initialise LDD simulation class
-simulation = ldd.LDDsimulation(tol=1E-14, debug=False, LDDsolver_tol=1E-9)
+simulation = ldd.LDDsimulation(tol=1E-14, debug=False, LDDsolver_tol=1E-6)
simulation.set_parameters(output_dir="./output/",
subdomain_def_points=subdomain_def_points,
isRichards=isRichards,
diff --git a/RR-multi-patch-with-inner-patch/RR-multi-patch-with-inner-patch.py b/RR-multi-patch-with-inner-patch/RR-multi-patch-with-inner-patch.py
index bbc1cb765265440a0c18f972702b5ad10bd7fcbb..f9d20f89bf9971227e3af55ad51516cedd415567 100755
--- a/RR-multi-patch-with-inner-patch/RR-multi-patch-with-inner-patch.py
+++ b/RR-multi-patch-with-inner-patch/RR-multi-patch-with-inner-patch.py
@@ -15,12 +15,12 @@ sym.init_printing()
# ----------------------------------------------------------------------------#
# ------------------- MESH ---------------------------------------------------#
# ----------------------------------------------------------------------------#
-mesh_resolution = 50
+mesh_resolution = 51
# ----------------------------------------:-------------------------------------#
# ------------------- TIME ---------------------------------------------------#
# ----------------------------------------------------------------------------#
timestep_size = 0.01
-number_of_timesteps = 500
+number_of_timesteps = 160
# 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
@@ -194,11 +194,11 @@ porosity = {
# subdom_num : subdomain L for L-scheme
L = {
- 1: {'wetting': 0.6},
- 2: {'wetting': 0.6},
- 3: {'wetting': 0.6},
- 4: {'wetting': 0.6},
- 5: {'wetting': 0.6}
+ 1: {'wetting': 0.25},
+ 2: {'wetting': 0.25},
+ 3: {'wetting': 0.25},
+ 4: {'wetting': 0.25},
+ 5: {'wetting': 0.25}
}
lamdal_w = 32
@@ -569,7 +569,7 @@ write_to_file = {
}
# initialise LDD simulation class
-simulation = ldd.LDDsimulation(tol=1E-14, debug=False, LDDsolver_tol=1E-9)
+simulation = ldd.LDDsimulation(tol=1E-14, debug=False, LDDsolver_tol=1E-7)
simulation.set_parameters(output_dir="./output/",
subdomain_def_points=subdomain_def_points,
isRichards=isRichards,
diff --git a/TP-TP-layered-soil-case-with-inner-patch-constant-solution/TP-TP-layered_soil_with_inner_patch_const_solution.py b/TP-TP-layered-soil-case-with-inner-patch-constant-solution/TP-TP-layered_soil_with_inner_patch_const_solution.py
index 4f469ec029df53d19985c6ae0ff63952a0c52f8b..8200e0e75521817d841efeee661a77b2572037d4 100755
--- a/TP-TP-layered-soil-case-with-inner-patch-constant-solution/TP-TP-layered_soil_with_inner_patch_const_solution.py
+++ b/TP-TP-layered-soil-case-with-inner-patch-constant-solution/TP-TP-layered_soil_with_inner_patch_const_solution.py
@@ -23,12 +23,12 @@ sym.init_printing()
# ----------------------------------------------------------------------------#
# ------------------- MESH ---------------------------------------------------#
# ----------------------------------------------------------------------------#
-mesh_resolution = 41
+mesh_resolution = 51
# ----------------------------------------:-----------------------------------#
# ------------------- TIME ---------------------------------------------------#
# ----------------------------------------------------------------------------#
-timestep_size = 0.00001
-number_of_timesteps = 100
+timestep_size = 0.01
+number_of_timesteps = 10
# 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 = 10
@@ -37,8 +37,8 @@ starttime = 0
Lw = 0.25
Lnw = 0.25
-l_param_w = 100
-l_param_nw = 100
+l_param_w = 200
+l_param_nw = 200
# global domain
subdomain0_vertices = [df.Point(-1.0,-1.0), #
@@ -510,103 +510,104 @@ ka_prime = {
}
+# # S-pc-relation ship. We use the van Genuchten approach, i.e. pc = 1/alpha*(S^{-1/m} -1)^1/n, where
+# # we set alpha = 0, assume m = 1-1/n (see Helmig) and assume that residual saturation is Sw
+# # this function needs to be monotonically decreasing in the capillary pressure pc.
+# # since in the richards case pc=-pw, this becomes as a function of pw a mono
+# # tonically INCREASING function like in our Richards-Richards paper. However
+# # since we unify the treatment in the code for Richards and two-phase, we need
+# # the same requierment
+# # for both cases, two-phase and Richards.
+# def saturation(pc, index, alpha):
+# # inverse capillary pressure-saturation-relationship
+# return df.conditional(pc > 0, 1/((1 + (alpha*pc)**index)**((index - 1)/index)), 1)
+# #
+# # S-pc-relation ship. We use the van Genuchten approach, i.e. pc = 1/alpha*(S^{-1/m} -1)^1/n, where
+# # we set alpha = 0, assume m = 1-1/n (see Helmig) and assume that residual saturation is Sw
+# def saturation_sym(pc, index, alpha):
+# # inverse capillary pressure-saturation-relationship
+# #df.conditional(pc > 0,
+# return 1/((1 + (alpha*pc)**index)**((index - 1)/index))
+#
+#
+# # derivative of S-pc relationship with respect to pc. This is needed for the
+# # construction of a analytic solution.
+# def saturation_sym_prime(pc, index, alpha):
+# # inverse capillary pressure-saturation-relationship
+# return -(alpha*(index - 1)*(alpha*pc)**(index - 1)) / ( (1 + (alpha*pc)**index)**((2*index - 1)/index) )
+#
+# # note that the conditional definition of S-pc in the nonsymbolic part will be
+# # incorporated in the construction of the exact solution below.
+# S_pc_sym = {
+# 1: ft.partial(saturation_sym, index=3, alpha=0.001),
+# 2: ft.partial(saturation_sym, index=3, alpha=0.001),
+# 3: ft.partial(saturation_sym, index=3, alpha=0.001),
+# 4: ft.partial(saturation_sym, index=3, alpha=0.001),
+# 5: ft.partial(saturation_sym, index=3, alpha=0.001),
+# 6: ft.partial(saturation_sym, index=3, alpha=0.001)
+# }
+#
+# S_pc_sym_prime = {
+# 1: ft.partial(saturation_sym_prime, index=3, alpha=0.001),
+# 2: ft.partial(saturation_sym_prime, index=3, alpha=0.001),
+# 3: ft.partial(saturation_sym_prime, index=3, alpha=0.001),
+# 4: ft.partial(saturation_sym_prime, index=3, alpha=0.001),
+# 5: ft.partial(saturation_sym_prime, index=3, alpha=0.001),
+# 6: ft.partial(saturation_sym_prime, index=3, alpha=0.001)
+# }
+#
+# sat_pressure_relationship = {
+# 1: ft.partial(saturation, index=3, alpha=0.001),
+# 2: ft.partial(saturation, index=3, alpha=0.001),
+# 3: ft.partial(saturation, index=3, alpha=0.001),
+# 4: ft.partial(saturation, index=3, alpha=0.001),
+# 5: ft.partial(saturation, index=3, alpha=0.001),
+# 6: ft.partial(saturation, index=3, alpha=0.001)
+# }
-# S-pc-relation ship. We use the van Genuchten approach, i.e. pc = 1/alpha*(S^{-1/m} -1)^1/n, where
-# we set alpha = 0, assume m = 1-1/n (see Helmig) and assume that residual saturation is Sw
-# this function needs to be monotonically decreasing in the capillary pressure pc.
-# since in the richards case pc=-pw, this becomes as a function of pw a mono
-# tonically INCREASING function like in our Richards-Richards paper. However
-# since we unify the treatment in the code for Richards and two-phase, we need
-# the same requierment
-# for both cases, two-phase and Richards.
-def saturation(pc, n_index, alpha):
+def saturation(pc, index):
# inverse capillary pressure-saturation-relationship
- return df.conditional(pc > 0, 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index)), 1)
-#
-# S-pc-relation ship. We use the van Genuchten approach, i.e. pc = 1/alpha*(S^{-1/m} -1)^1/n, where
-# we set alpha = 0, assume m = 1-1/n (see Helmig) and assume that residual saturation is Sw
-def saturation_sym(pc, n_index, alpha):
+ return df.conditional(pc > 0, 1/((1 + pc)**(1/(index + 1))), 1)
+
+
+def saturation_sym(pc, index):
# inverse capillary pressure-saturation-relationship
- #df.conditional(pc > 0,
- return 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index))
+ return 1/((1 + pc)**(1/(index + 1)))
# derivative of S-pc relationship with respect to pc. This is needed for the
# construction of a analytic solution.
-def saturation_sym_prime(pc, n_index, alpha):
+def saturation_sym_prime(pc, index):
# inverse capillary pressure-saturation-relationship
- return -(alpha*(n_index - 1)*(alpha*pc)**(n_index - 1)) / ( (1 + (alpha*pc)**n_index)**((2*n_index - 1)/n_index) )
-
-# note that the conditional definition of S-pc in the nonsymbolic part will be
-# incorporated in the construction of the exact solution below.
+ return -1/((index+1)*(1 + pc)**((index+2)/(index+1)))
+#
S_pc_sym = {
- 1: ft.partial(saturation_sym, n_index=3, alpha=0.001),
- 2: ft.partial(saturation_sym, n_index=3, alpha=0.001),
- 3: ft.partial(saturation_sym, n_index=3, alpha=0.001),
- 4: ft.partial(saturation_sym, n_index=3, alpha=0.001),
- 5: ft.partial(saturation_sym, n_index=3, alpha=0.001),
- 6: ft.partial(saturation_sym, n_index=3, alpha=0.001)
+ 1: ft.partial(saturation_sym, index=1),
+ 2: ft.partial(saturation_sym, index=1),
+ 3: ft.partial(saturation_sym, index=1),
+ 4: ft.partial(saturation_sym, index=1),
+ 5: ft.partial(saturation_sym, index=1),
+ 6: ft.partial(saturation_sym, index=1)
}
S_pc_sym_prime = {
- 1: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
- 2: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
- 3: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
- 4: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
- 5: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
- 6: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001)
+ 1: ft.partial(saturation_sym_prime, index=1),
+ 2: ft.partial(saturation_sym_prime, index=1),
+ 3: ft.partial(saturation_sym_prime, index=1),
+ 4: ft.partial(saturation_sym_prime, index=1),
+ 5: ft.partial(saturation_sym_prime, index=1),
+ 6: ft.partial(saturation_sym_prime, index=1)
}
sat_pressure_relationship = {
- 1: ft.partial(saturation, n_index=3, alpha=0.001),
- 2: ft.partial(saturation, n_index=3, alpha=0.001),
- 3: ft.partial(saturation, n_index=3, alpha=0.001),
- 4: ft.partial(saturation, n_index=3, alpha=0.001),
- 5: ft.partial(saturation, n_index=3, alpha=0.001),
- 6: ft.partial(saturation, n_index=3, alpha=0.001)
+ 1: ft.partial(saturation, index=1),
+ 2: ft.partial(saturation, index=1),
+ 3: ft.partial(saturation, index=1),
+ 4: ft.partial(saturation, index=1),
+ 5: ft.partial(saturation, index=1),
+ 6: ft.partial(saturation, index=1)
}
-# def saturation(pc, n_index):
-# # inverse capillary pressure-saturation-relationship
-# return df.conditional(pc > 0, 1/((1 + pc)**(1/(n_index + 1))), 1)
-#
-#
-# def saturation_sym(pc, n_index):
-# # inverse capillary pressure-saturation-relationship
-# return 1/((1 + pc)**(1/(n_index + 1)))
-#
-# def saturation_sym_prime(pc, n_index):
-# # inverse capillary pressure-saturation-relationship
-# return -1/((n_index+1)*(1 + pc)**((n_index+2)/(n_index+1)))
-#
-#
-# S_pc_sym = {
-# 1: ft.partial(saturation_sym, n_index=1),
-# 2: ft.partial(saturation_sym, n_index=1),
-# 3: ft.partial(saturation_sym, n_index=1),
-# 4: ft.partial(saturation_sym, n_index=1),
-# 5: ft.partial(saturation_sym, n_index=1),
-# 6: ft.partial(saturation_sym, n_index=1)
-# }
-#
-# S_pc_sym_prime = {
-# 1: ft.partial(saturation_sym_prime, n_index=1),
-# 2: ft.partial(saturation_sym_prime, n_index=1),
-# 3: ft.partial(saturation_sym_prime, n_index=1),
-# 4: ft.partial(saturation_sym_prime, n_index=1),
-# 5: ft.partial(saturation_sym_prime, n_index=1),
-# 6: ft.partial(saturation_sym_prime, n_index=1)
-# }
-#
-# sat_pressure_relationship = {
-# 1: ft.partial(saturation, n_index=1),
-# 2: ft.partial(saturation, n_index=1),
-# 3: ft.partial(saturation, n_index=1),
-# 4: ft.partial(saturation, n_index=1),
-# 5: ft.partial(saturation, n_index=1),
-# 6: ft.partial(saturation, n_index=1)
-# }
-
#############################################
# Manufacture source expressions with sympy #
@@ -627,7 +628,7 @@ for subdomain, isR in isRichards.items():
p_e_sym[subdomain].update({phase: p_dict[phase]})
if isR:
- pc_e_sym.update({subdomain: -p_e_sym[subdomain]['wetting']})
+ pc_e_sym.update({subdomain: -p_e_sym[subdomain]['wetting'].copy()})
else:
pc_e_sym.update({subdomain: p_e_sym[subdomain]['nonwetting'] - p_e_sym[subdomain]['wetting']})