diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index 4699e0e493d6450c73ebca7f2515b2fcffbb89f5..41c91bb0e4597f5c41fcf959e462a2ed1bfb88ae 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -382,7 +382,10 @@ class LDDsimulation(object):
analyse_condition=False)
# calculate spacetime errornorms for self.output.
- space_errornorm = self._spacetime_errornorm_calc_step()
+ if self.exact_solution is not None:
+ space_errornorm = self._spacetime_errornorm_calc_step()
+ else:
+ space_errornorm = None
# write solutions to files. It is paramount, that self.timestep_num
# be only updated after calling self.write_data_to_disc()
@@ -535,13 +538,6 @@ class LDDsimulation(object):
subdomain.pressure_prev_iter[phase].assign(
subdomain.pressure[phase]
)
- # onyl becomes active in TPR case
- if self.include_gravity and subdomain.isRichards and subdomain.TPR_occures:
- for interface in subdomain.has_interface_with_TP_neighbour:
- subdomain.calc_gravity_neumann_flux(
- interface_index=interface,
- phase="nonwetting"
- )
# end loop over subdomains
##### stopping criterion for the solver.
total_subsequent_error = np.amax(max_subsequent_error)
@@ -652,17 +648,6 @@ class LDDsimulation(object):
subdomain.write_pressure_to_interfaces()
subdomain.write_pressure_to_interfaces(previous_iter=True)
subdomain.calc_gl0_term()
- # If gravity is included and a Richards subdomain has a TPR coupling,
- # we need to calculate the Neumann flux of the nonwetting phase
- # that appears due to gravity
- if self.include_gravity:
- if subdomain.isRichards and subdomain.TPR_occures:
- for interface in subdomain.has_interface_with_TP_neighbour:
- subdomain.calc_gravity_neumann_flux(
- initial=True,
- interface_index=interface,
- phase="nonwetting"
- )
return solution_over_iteration_within_timestep
diff --git a/LDDsimulation/boundary_and_interface.py b/LDDsimulation/boundary_and_interface.py
index 7d1a962d09fe354abd30d4814108479472acd430..121515bd506b46046b517e9eecfcb61b63eec4ce 100644
--- a/LDDsimulation/boundary_and_interface.py
+++ b/LDDsimulation/boundary_and_interface.py
@@ -615,8 +615,6 @@ class Interface(BoundaryPart):
else: #function == "gli":
self.gli_term.update({subdomain_index: dict()})
self.gli_term_prev.update({subdomain_index: dict()})
- self.gravity_neumann_flux.update({subdomain_index: dict()})
- self.gravity_neumann_flux_prev.update({subdomain_index: dict()})
for phase in has_phases:
if function == "pressure":
@@ -625,8 +623,7 @@ class Interface(BoundaryPart):
else: #function == "gli":
self.gli_term[subdomain_index].update({phase: dict()})
self.gli_term_prev[subdomain_index].update({phase: dict()})
- self.gravity_neumann_flux[subdomain_index].update({phase: dict()})
- self.gravity_neumann_flux_prev[subdomain_index].update({phase: dict()})
+
space = function_space[function][phase]
dof_coordinates = self.dofs_and_coordinates(space, interface_marker, interface_marker_value)
# the dof dictionaries are grouped by the global (with respect to the
@@ -640,8 +637,6 @@ class Interface(BoundaryPart):
else: #function == "gli":
self.gli_term[subdomain_index][phase].update({global_mesh_facet_index: dict()})
self.gli_term_prev[subdomain_index][phase].update({global_mesh_facet_index: dict()})
- self.gravity_neumann_flux[subdomain_index][phase].update({global_mesh_facet_index: dict()})
- self.gravity_neumann_flux_prev[subdomain_index][phase].update({global_mesh_facet_index: dict()})
for dof_coord in facet_dof_coord_dict.values():
# we want to use the coordinates as key in a dictionary.
@@ -655,8 +650,6 @@ class Interface(BoundaryPart):
else: #function == "gli":
self.gli_term[subdomain_index][phase][global_mesh_facet_index].update({coord_tupel: 0.0})
self.gli_term_prev[subdomain_index][phase][global_mesh_facet_index].update({coord_tupel: 0.0})
- self.gravity_neumann_flux[subdomain_index][phase][global_mesh_facet_index].update({coord_tupel: 0.0})
- self.gravity_neumann_flux_prev[subdomain_index][phase][global_mesh_facet_index].update({coord_tupel: 0.0})
# end init_interface_values
diff --git a/LDDsimulation/domainPatch.py b/LDDsimulation/domainPatch.py
index da2596a07b066e6f6e0c8809e21d91ce4473d73d..248b43216995106bae8194de0ccb6770543a1c4e 100644
--- a/LDDsimulation/domainPatch.py
+++ b/LDDsimulation/domainPatch.py
@@ -248,14 +248,6 @@ class DomainPatch(df.SubDomain):
self._calc_interface_dof_indices_and_coordinates()
self._calc_corresponding_dof_indices()
self._init_interface_values()
- # # Dictionary of FEM function holding the same dofs gli_assembled above
- # # but as a df.Function. This is needed for writing out the gli terms when
- # # analyising the convergence of a timestep.
- # self.gli_function = dict()
- # for phase in self.has_phases:
- # self.gli_function.update(
- # {phase: df.Function(self.function_space["gli"][phase])}
- # )
# END constructor
@@ -373,16 +365,6 @@ class DomainPatch(df.SubDomain):
# gather interface forms
interface_forms = []
gli_forms = []
- # the following needs to be executed on a TP domain with R neighbours
- # in case gravity is included.
- if self.include_gravity and phase=="nonwetting" and self.TPR_occures:
- rhs_gravity_term_interface = []
- for interface in self.has_interface_with_Richards_neighbour:
- marker_value = self.interface[interface].marker_value
- gravity_neumann_flux = self.read_gravity_neumann_flux(interface_index=interface, phase="nonwetting")
- rhs_gravity_term_interface.append(
- -dt*gravity_neumann_flux*phi_a*ds(marker_value)
- )
# 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")
@@ -427,12 +409,8 @@ class DomainPatch(df.SubDomain):
if self.has_interface is not None:
form = form1 + form2 + sum(interface_forms)
if self.include_gravity:
- if phase=="nonwetting" and self.TPR_occures:
- rhs = rhs1 + rhs2 + dt*source_a*phi_a*dx - sum(gli_forms) \
- + rhs_gravity + sum(rhs_gravity_term_interface)
- else:
- # z_a included the minus sign already, see self._calc_gravity_expressions
- rhs = rhs1 + rhs2 + dt*source_a*phi_a*dx - sum(gli_forms) + rhs_gravity
+ # z_a included the minus sign already, see self._calc_gravity_expressions
+ rhs = rhs1 + rhs2 + dt*source_a*phi_a*dx - sum(gli_forms) + rhs_gravity
else:
rhs = rhs1 + rhs2 + dt*source_a*phi_a*dx - sum(gli_forms)
else:
@@ -449,198 +427,6 @@ class DomainPatch(df.SubDomain):
return form_and_rhs
- def calc_gravity_neumann_flux(self,
- interface_index: int,
- phase: str,
- initial: bool = False,
- debug: bool = False
- ) -> None:
- """calculate the gravity term of the neumann flux
-
- This method is written such that it can calculate the gravity neumann
- flux for both phases in case it is needed, but is actually only used
- for the flux of the nonwetting gravity term of a Richards domain with
- TPR coupling.
- """
- subdomain = self.subdomain_index
- dt = self.timestep_size
- Lambda = self.lambda_param
- interf_ind = interface_index
-
- interface = self.interface[interf_ind]
-
- if initial:
- pressure_previous = self.pressure_prev_timestep
- else:
- pressure_previous = self.pressure_prev_iter
-
- if self.isRichards:
- # assuming that we normalised atmospheric pressure to zero,
- # pc = -pw in the Richards case.
- pc_prev_iter = -pressure_previous['wetting']
- else:
- pw_prev_iter = pressure_previous['wetting']
- pnw_prev_iter = pressure_previous['nonwetting']
- pc_prev_iter = pnw_prev - pw_prev
-
- # n = df.FacetNormal(self.mesh)
- S = self.saturation
-
- # viscosity
- mu_a = self.viscosity[phase]
- # relative permeability
- ka = self.relative_permeability[phase]
- z_a = self.gravity_term[phase]
-
- if phase == 'wetting':
- # z_a included the minus sign, see self._calc_gravity_expressions
- flux = 1/mu_a*ka(S(pc_prev_iter))*df.grad(z_a)
- else:
- # phase == 'nonwetting'
- flux = 1/mu_a*ka(1-S(pc_prev_iter))*df.grad(z_a)
-
- flux_function = df.project(flux, self.function_space["flux"][phase])
- flux_x, flux_y = flux_function.split()
- gravity_neumann_flux = 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.
- neumann_flux_dof_coords = 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():
- neumann_flux_dofs_along_facet = neumann_flux_dof_coords[local_facet_ind]
- global_facet_ind = local2global_facet[local_facet_ind]
- for neumann_flux_dof_index, neumann_flux_dof_coord in neumann_flux_dofs_along_facet.items():
- flux_x_dof_index = corresponding_dof_index[local_facet_ind][neumann_flux_dof_index]["flux_x"]
- flux_y_dof_index = corresponding_dof_index[local_facet_ind][neumann_flux_dof_index]["flux_y"]
-
- gravity_neumann_flux.vector()[neumann_flux_dof_index] = flux_x.vector()[flux_x_dof_index]*normal[0] \
- + flux_y.vector()[flux_y_dof_index]*normal[1]
-
- # save this newly calculated dof to the interface
- # dictionary neumann_flux_save_dict
- dof_coord_tupel = (neumann_flux_dof_coord[0], neumann_flux_dof_coord[1])
- if initial:
- # because we don't know if the Richards domain or the TP domain
- # will be up next in the solver, we need to save the initial
- # gravity neumann flux to both the current and the prev dictionary.
- interface.gravity_neumann_flux_prev[subdomain][phase][global_facet_ind].update(
- {dof_coord_tupel: gravity_neumann_flux.vector()[neumann_flux_dof_index]}
- )
- interface.gravity_neumann_flux[subdomain][phase][global_facet_ind].update(
- {dof_coord_tupel: gravity_neumann_flux.vector()[neumann_flux_dof_index]}
- )
- else:
- # in the noninitial case we always save to the current dictionary
- # similarly to the pressures.
- interface.gravity_neumann_flux[subdomain][phase][global_facet_ind].update(
- {dof_coord_tupel: gravity_neumann_flux.vector()[neumann_flux_dof_index]}
- )
-
-
- def read_gravity_neumann_flux(self,
- interface_index: int,
- phase: str,
- debug: bool = False) -> df.Function: # tp.Dict[str, fl.Form]
- """in the case of a TPR coupling read the additional gravity neumann flux
- for the nonwetting phase from the neighbour
-
- IMPLEMENTATION
- same logic as in the calc_gli_term method applies.
-
- PARAMETERS
- ------------------------------------
- interface_index int interface index for which to calculate the
- gravity_neumann_flux term.
- phase str phase to calculate the gravity_neumann_flux term for.
- # DEBUG: bool toggle the debug output.
- """
- # this is the number of the current iteration of the subdomain.
- iteration = self.iteration_number
- subdomain = self.subdomain_index
- interface = self.interface[interface_index]
- # gravity_neumann_flux should be initialised by zero.
- gravity_neumann_flux = df.Function(self.function_space["gli"][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]
-
- # read previous neighbouring gravity neumann flux 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.gravity_neumann_flux[neighbour] holds the "previous neumann gravity flux"
- # we need to read these values from the interface.
- gravity_flux_read_dict = interface.gravity_neumann_flux[neighbour][phase]
-
- else:
- # 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.gravity_neumann_flux_prev holds the "previous neumann gravity flux"
- # we need to read from the interface.
- gravity_flux_read_dict = interface.gravity_neumann_flux_prev[neighbour][phase]
-
- flux_interface_dofs = self.interface_dof_indices_and_coordinates["gli"][interface_index][phase]
- local2global_facet = interface.local_to_global_facet_map[subdomain]
- for local_facet_index, flux_facet_dofs2coordinates in flux_interface_dofs.items():
- global_facet_index = local2global_facet[local_facet_index]
- # read gravity_neumann_flux_prev term from the neighbour.
- neighbour_flux_on_facet = gravity_flux_read_dict[global_facet_index]
-
- for flux_dof_index, flux_dof_coord in flux_facet_dofs2coordinates.items():
- # same with the gravity_neumann_flux_previous_dof
- gravity_neumann_flux_neighbour_dof = None
- for flux_neighbour_coord_tupel, flux_neighbour_dof in neighbour_flux_on_facet.items():
- flux_neighbour_coord = np.array([flux_neighbour_coord_tupel[0], flux_neighbour_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(flux_dof_coord, flux_neighbour_coord, rtol=1e-10, atol=1e-14)
- b_close_a = np.allclose(flux_neighbour_coord, flux_dof_coord, rtol=1e-10, atol=1e-14)
- if a_close_b and b_close_a:
- gravity_neumann_flux_neighbour_dof = flux_neighbour_dof
- break
-
- if gravity_neumann_flux_neighbour_dof is None:
- raise RuntimeError(f"didn't find gravity_neumann_flux_neighbour_dof corresponding to dict key ({flux_dof_coord})",
- "something is wrong")
-
- gravity_neumann_flux.vector()[flux_dof_index] = gravity_neumann_flux_neighbour_dof
-
- # In the following case, the gravity_neumann_flux of the neighbour is saved to
- # gravity_neumann_flux currently and will be overwritten in the next step,
- # if we don't save it to gravity_neumann_flux_prev of the neighbour.
- 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.gravity_neumann_flux_prev[neighbour][phase].update(#
- {global_facet_index: interface.gravity_neumann_flux[neighbour][phase][global_facet_index].copy()}
- )
-
- return gravity_neumann_flux
- ### END read_gravity_neumann_flux
-
-
def calc_gl0_term(self, debug: bool = False) -> None: # tp.Dict[str, fl.Form]
"""calculate the gl0 terms for all interfaces of the subdomain.
@@ -671,10 +457,10 @@ class DomainPatch(df.SubDomain):
S = self.saturation
phases_for_gl0 = self.has_phases
- # if self.isRichards and self.TPR_occures:
- # interface_has_TP_neighbour = not interface.neighbour_isRichards[self.subdomain_index]
- # if interface_has_TP_neighbour and self.include_gravity:
- # phases_for_gl0 = ['wetting', 'nonwetting']
+ if self.isRichards and self.TPR_occures:
+ interface_has_TP_neighbour = not interface.neighbour_isRichards[self.subdomain_index]
+ if interface_has_TP_neighbour and self.include_gravity:
+ phases_for_gl0 = ['wetting', 'nonwetting']
for phase in phases_for_gl0:
@@ -683,17 +469,18 @@ class DomainPatch(df.SubDomain):
# relative permeability
ka = self.relative_permeability[phase]
# previous timestep pressure
- # if self.isRichards and interface_has_TP_neighbour and phase is "nonwetting":
- # pa_prev = 0
- # else:
- pa_prev = self.pressure_prev_timestep[phase]
+ if self.isRichards and interface_has_TP_neighbour and phase is "nonwetting":
+ # dummy definition. we don't use this in this case
+ pa_prev = 0
+ else:
+ pa_prev = self.pressure_prev_timestep[phase]
if self.include_gravity:
z_a = self.gravity_term[phase]
+
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)
@@ -701,17 +488,16 @@ class DomainPatch(df.SubDomain):
flux = -1/mu_a*ka(S(pc_prev))*df.grad(pa_prev)
else:
# phase == 'nonwetting'
- # if self.isRichards and interface_has_TP_neighbour and phase is "nonwetting":
- # if self.include_gravity:
- # flux = -1/mu_a*ka(1-S(pc_prev))*df.grad(-z_a)
- # else:
- # flux = 0
- # else:
- # 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)
+ if self.isRichards and interface_has_TP_neighbour and phase is "nonwetting":
+ if self.include_gravity:
+ flux = -1/mu_a*ka(1-S(pc_prev))*df.grad(-z_a)
+ else:
+ flux = 0
else:
- flux = -1/mu_a*ka(1-S(pc_prev))*df.grad(pa_prev)
+ if self.include_gravity:
+ flux = -1/mu_a*ka(1-S(pc_prev))*df.grad(pa_prev - z_a)
+ else:
+ flux = -1/mu_a*ka(1-S(pc_prev))*df.grad(pa_prev)
flux_function = df.project(flux, self.function_space["flux"][phase])
flux_x, flux_y = flux_function.split()
@@ -738,15 +524,15 @@ class DomainPatch(df.SubDomain):
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"]
-
- # if self.isRichards and interface_has_TP_neighbour and phase is "nonwetting" and self.include_gravity:
- # gl0.vector()[gli_dof_index] = flux_x.vector()[flux_x_dof_index]*normal[0] \
- # + flux_y.vector()[flux_y_dof_index]*normal[1]
- # else:
- 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 self.isRichards and interface_has_TP_neighbour and phase is "nonwetting" and self.include_gravity:
+ gl0.vector()[gli_dof_index] = flux_x.vector()[flux_x_dof_index]*normal[0] \
+ + flux_y.vector()[flux_y_dof_index]*normal[1]
+ else:
+ 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]}")
@@ -1356,12 +1142,12 @@ class DomainPatch(df.SubDomain):
marker_value = self.interface[interf].marker_value
space = self.function_space
- neighbour_assumes_R = self.interface[interf].neighbour_isRichards[self.subdomain_index]
+ neighbour_assumes_TP = not 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"]
+ if self.isRichards and neighbour_assumes_TP:
+ subdom_has_phases=["wetting", "nonwetting"]
else:
subdom_has_phases=self.has_phases
diff --git a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-pure-dd-realistic.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-pure-dd-realistic.py
old mode 100644
new mode 100755
index 3fe92a0bf76129aa0887b93ebf5c53c5a6c15d05..13d24f67a13880cf08c2c74295fe8cbaa3661830
--- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-pure-dd-realistic.py
+++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-pure-dd-realistic.py
@@ -30,27 +30,26 @@ resolutions = {
# 4: 7e-7,
# 8: 7e-7,
# 16: 7e-7,
- 32: 1e-6,
+ 32: 5e-6,
# 64: 7e-7,
# 128: 7e-7,
# 256: 7e-7
}
############ GRID #######################
-# mesh_resolution = 20
-timestep_size = 0.0005
-number_of_timesteps = 20
+timestep_size = 0.001
+number_of_timesteps = 15
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 = 1
+number_of_timesteps_to_analyse = 5
starttimes = [0.0]
-Lw = 10 #/timestep_size
-Lnw= 50
+Lw = 0.04 #/timestep_size
+Lnw= 0.8
-lambda_w = 40000
-lambda_nw = 40000
+lambda_w = 2
+lambda_nw = 8
include_gravity = True
debugflag = True
@@ -190,14 +189,21 @@ lambda_param = {#
'nonwetting': lambda_nw}
}
+intrinsic_permeability = {
+ 1: {"wetting": 10e-6,
+ "nonwetting": 10e-6},
+ 2: {"wetting": 10e-6,
+ "nonwetting": 10e-6},
+}
+
## relative permeabilty functions on subdomain 1
def rel_perm1w(s):
# relative permeabilty wetting on subdomain1
- return s**2
+ return intrinsic_permeability[1]["wetting"]*s**2
def rel_perm1nw(s):
# relative permeabilty nonwetting on subdomain1
- return (1-s)**2
+ return intrinsic_permeability[1]["nonwetting"]*(1-s)**2
_rel_perm1w = ft.partial(rel_perm1w)
_rel_perm1nw = ft.partial(rel_perm1nw)
@@ -232,11 +238,11 @@ relative_permeability = {#
# relative permeabilty functions on subdomain 1
def rel_perm1w_prime(s):
# relative permeabilty on subdomain1
- return 2*s
+ return intrinsic_permeability[1]["wetting"]*2*s
def rel_perm1nw_prime(s):
# relative permeabilty on subdomain1
- return -2*(1-s)
+ return -1*intrinsic_permeability[1]["nonwetting"]*2*(1-s)
# # definition of the derivatives of the relative permeabilities
# # relative permeabilty functions on subdomain 1
@@ -340,22 +346,22 @@ def saturation_sym_prime(pc, 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=2),
- 2: ft.partial(saturation_sym, index=2),
+ 1: ft.partial(saturation_sym, index=1),
+ 2: ft.partial(saturation_sym, index=1),
# 3: ft.partial(saturation_sym, index=2),
# 4: ft.partial(saturation_sym, index=1)
}
S_pc_sym_prime = {
- 1: ft.partial(saturation_sym_prime, index=2),
- 2: ft.partial(saturation_sym_prime, index=2),
+ 1: ft.partial(saturation_sym_prime, index=1),
+ 2: ft.partial(saturation_sym_prime, index=1),
# 3: ft.partial(saturation_sym_prime, index=2),
# 4: ft.partial(saturation_sym_prime, index=1)
}
sat_pressure_relationship = {
- 1: ft.partial(saturation, index=2),
- 2: ft.partial(saturation, index=2),
+ 1: ft.partial(saturation, index=1),
+ 2: ft.partial(saturation, index=1),
# 3: ft.partial(saturation, index=2),
# 4: ft.partial(saturation, index=1)
}
@@ -373,6 +379,7 @@ t = sym.symbols('t', positive=True)
# '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)
+
p_e_sym = {
1: {'wetting': (-6 - (1+t*t)*(1 + x*x + y*y)), #*cutoff,
'nonwetting': (-1 -t*(1.1+ y*y))}, #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2},
diff --git a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-pure-dd.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-pure-dd.py
index bd55ce3c7bc9f6a76d0cec9b3ffcc92402a937bc..2c861d2f327a3563950cefe4fc27a93e4dfbfdb8 100755
--- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-pure-dd.py
+++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-pure-dd.py
@@ -30,30 +30,30 @@ resolutions = {
# 4: 7e-7,
# 8: 7e-7,
# 16: 7e-7,
- 32: 1e-6,
- # 64: 7e-7,
+ # 32: 1e-6,
+ 64: 7e-7,
# 128: 7e-7,
# 256: 7e-7
}
############ GRID #######################
# mesh_resolution = 20
-timestep_size = 0.01
-number_of_timesteps = 150
+timestep_size = 0.005
+number_of_timesteps = 250
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
starttimes = [0.0]
-Lw = 0.1 #/timestep_size
-Lnw= 0.1
+Lw = 0.025 #/timestep_size
+Lnw= 0.025
-lambda_w = 4
-lambda_nw = 4
+lambda_w = 2
+lambda_nw = 2
include_gravity = True
debugflag = False
-analyse_condition = False
+analyse_condition = True
output_string = "./output/{}-{}_timesteps{}_P{}".format(datestr, use_case, number_of_timesteps, FEM_Lagrange_degree)
@@ -170,7 +170,7 @@ densities = {
'nonwetting': 1} #1.225},
}
-gravity_acceleration = 9.81
+gravity_acceleration = 1 #9.81
L = {#
# subdom_num : subdomain L for L-scheme
@@ -373,10 +373,10 @@ t = sym.symbols('t', positive=True)
# } #-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)
p_e_sym = {
- 1: {'wetting': (-6 - (1+t*t)*(1 + x*x + y*y)), #*cutoff,
- 'nonwetting': (-1 -t*(1.1+ y*y))}, #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2},
- 2: {'wetting': (-6.0 - (1.0 + t*t)*(1.0 + x*x)), #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2,
- 'nonwetting': (-1 -t*(1.1 + y*y) - sym.sin((x*y-0.5*t)*y**2)**2)}, #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2},
+ 1: {'wetting': (-6 - (1+t*t)*(1 + x*x + y*y))}, #*cutoff,
+ # 'nonwetting': (-1 -t*(1.1+ y*y))}, #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2},
+ 2: {'wetting': (-6 - (1+t*t)*(1 + x*x + y*y)), #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2,
+ 'nonwetting': -t*(1.1 + y*x)*(sym.sin((x*y-0.5*t)*y**3))**4}, #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2},
}
pc_e_sym = dict()
diff --git a/Two-phase-Two-phase/two-patch/injection/TP-TP-2-patch-injection.py b/Two-phase-Two-phase/two-patch/injection/TP-TP-2-patch-injection.py
index 6206cd934a9cc9d2930cc7c73069686509c99fe4..ef3b3d573fa562dd45170a17c9f22b4a131a69c1 100755
--- a/Two-phase-Two-phase/two-patch/injection/TP-TP-2-patch-injection.py
+++ b/Two-phase-Two-phase/two-patch/injection/TP-TP-2-patch-injection.py
@@ -38,7 +38,7 @@ resolutions = {
############ GRID #######################
# mesh_resolution = 20
-timestep_size = 0.000001
+timestep_size = 0.0005
number_of_timesteps = 13
plot_timestep_every = 1
# decide how many timesteps you want analysed. Analysed means, that we write out
@@ -46,11 +46,11 @@ plot_timestep_every = 1
number_of_timesteps_to_analyse = 0
starttime = 0.0
-Lw = 0.5 #/timestep_size
+Lw = 0.05 #/timestep_size
Lnw=Lw
-lambda_w = 4000
-lambda_nw = 4000
+lambda_w = 80
+lambda_nw = 80
include_gravity = False
debugflag = True
@@ -177,6 +177,14 @@ densities = {
'nonwetting': 1225}, #1225},
}
+intrinsic_permeability = {
+ 1: {"wetting": 10e-6,
+ "nonwetting": 10e-6},
+ 2: {"wetting": 10e-6,
+ "nonwetting": 10e-6},
+}
+
+
gravity_acceleration = 9.81
L = {#
@@ -199,11 +207,11 @@ lambda_param = {#
## relative permeabilty functions on subdomain 1
def rel_perm1w(s):
# relative permeabilty wetting on subdomain1
- return s**2
+ return intrinsic_permeability[1]["wetting"]*s**2
def rel_perm1nw(s):
# relative permeabilty nonwetting on subdomain1
- return (1-s)**2
+ return intrinsic_permeability[1]["nonwetting"]*(1-s)**2
_rel_perm1w = ft.partial(rel_perm1w)
_rel_perm1nw = ft.partial(rel_perm1nw)
@@ -215,10 +223,10 @@ subdomain1_rel_perm = {
## relative permeabilty functions on subdomain 2
def rel_perm2w(s):
# relative permeabilty wetting on subdomain2
- return s**3
+ return intrinsic_permeability[2]["wetting"]*s**3
def rel_perm2nw(s):
# relative permeabilty nonwetting on subdosym.cos(0.8*t - (0.8*x + 1/7*y))main2
- return (1-s)**3
+ return intrinsic_permeability[2]["nonwetting"]*(1-s)**3
_rel_perm2w = ft.partial(rel_perm2w)
_rel_perm2nw = ft.partial(rel_perm2nw)
@@ -239,21 +247,21 @@ relative_permeability = {#
# relative permeabilty functions on subdomain 1
def rel_perm1w_prime(s):
# relative permeabilty on subdomain1
- return 2*s
+ return intrinsic_permeability[1]["wetting"]*2*s
def rel_perm1nw_prime(s):
# relative permeabilty on subdomain1
- return -2*(1-s)
+ return -1*intrinsic_permeability[1]["nonwetting"]*2*(1-s)
# # definition of the derivatives of the relative permeabilities
# # relative permeabilty functions on subdomain 1
def rel_perm2w_prime(s):
# relative permeabilty on subdomain1
- return 3*s**2
+ return intrinsic_permeability[2]["wetting"]*3*s**2
def rel_perm2nw_prime(s):
# relative permeabilty on subdomain1
- return -3*(1-s)**2
+ return -3*intrinsic_permeability[2]["nonwetting"]*(1-s)**2
_rel_perm1w_prime = ft.partial(rel_perm1w_prime)
_rel_perm1nw_prime = ft.partial(rel_perm1nw_prime)
@@ -378,8 +386,8 @@ t = sym.symbols('t', positive=True)
initial_condition = {
1: {'wetting': sym.printing.ccode(-6 - (1+t*t)*(1 + x*x + y*y)), #*cutoff,
'nonwetting': sym.printing.ccode(-1 -t*(1.1+ y*y))}, #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2},
- 2: {'wetting': sym.printing.ccode(-6.0 - (1.0 + t*t)*(1.0 + x*x)), #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2,
- 'nonwetting': sym.printing.ccode(-1 -t*(1.1 + y*y) - sym.sin((x*y-0.5*t)*y**2)**2)}, #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2},
+ 2: {'wetting': sym.printing.ccode(-6 - (1+t*t)*(1 + x*x + y*y)), #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2,
+ 'nonwetting': sym.printing.ccode(-1 -t*(1.1 + y*y))}, #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2},
}
### constructing source experessions.
@@ -408,7 +416,7 @@ extraction_radius = 0.075
#
def mollifier2d(x, y, epsilon):
""" one d mollifier """
- out_expr = sym.exp(-1/(1-(x**2 + y**2)/epsilon**2))
+ out_expr = 0.05*sym.exp(-1/(1-(x**2 + y**2)/epsilon**2))
return out_expr
mollifier2d_handle_i = ft.partial(mollifier2d, epsilon=injection_radius)