diff --git a/LDDsimulation/domainPatch.py b/LDDsimulation/domainPatch.py
index 8bb52d4ec05f626ec47efe98ed74a4b80f2ff3c8..099fc74e82c08b7a4ca79330bdd9ee16102ec27f 100644
--- a/LDDsimulation/domainPatch.py
+++ b/LDDsimulation/domainPatch.py
@@ -176,7 +176,7 @@ class DomainPatch(df.SubDomain):
self.has_phases = ['wetting']
else:
self.has_phases = ['wetting', 'nonwetting'] #['nonwetting', 'wetting']
-
+
if self.has_interface is not None:
if self.isRichards:
self.has_interface_with_Richards_neighbour = None
@@ -491,7 +491,7 @@ 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.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:
diff --git a/LDDsimulation/helpers.py b/LDDsimulation/helpers.py
index 5cf02e8fe3f038bfa2804e7160ae6373bd648e2b..c73f81b67ade74348adc30f838ae6d4e860ef80a 100644
--- a/LDDsimulation/helpers.py
+++ b/LDDsimulation/helpers.py
@@ -111,23 +111,23 @@ def generate_exact_solution_expressions(
if phase == "nonwetting":
# x part of div(flux) for nonwetting
- dxdxflux = -1/mu*dka(1-S)*dS*dxpc*dxpa + 1/mu*dxdxpa*ka(1-S)
+ div_flux_x = -1*(-1/mu*dka(1-S)*dS*dxpc*dxpa + 1/mu*dxdxpa*ka(1-S))
# y part of div(flux) for nonwetting
if include_gravity:
- dydyflux = -1/mu*dka(1-S)*dS*dypc*(dypa + rho*g) + 1/mu*dydypa*ka(1-S)
+ div_flux_y = -1*(-1/mu*dka(1-S)*dS*dypc*(dypa + rho*g) + 1/mu*dydypa*ka(1-S))
else:
- dydyflux = -1/mu*dka(1-S)*dS*dypc*dypa + 1/mu*dydypa*ka(1-S)
+ div_flux_y = -1*(-1/mu*dka(1-S)*dS*dypc*dypa + 1/mu*dydypa*ka(1-S))
else:
# x part of div(flux) for wetting
- dxdxflux = 1/mu*dka(S)*dS*dxpc*dxpa + 1/mu*dxdxpa*ka(S)
+ div_flux_x = -1*(1/mu*dka(S)*dS*dxpc*dxpa + 1/mu*dxdxpa*ka(S))
# y part of div(flux) for wetting
if include_gravity:
- dydyflux = 1/mu*dka(S)*dS*dypc*(dypa + rho*g) + 1/mu*dydypa*ka(S)
+ div_flux_y = -1*(1/mu*dka(S)*dS*dypc*(dypa + rho*g) + 1/mu*dydypa*ka(S))
else:
- dydyflux = 1/mu*dka(S)*dS*dypc*(dypa) + 1/mu*dydypa*ka(S)
+ div_flux_y = -1*(1/mu*dka(S)*dS*dypc*(dypa) + 1/mu*dydypa*ka(S))
- div_flux[subdomain].update({phase: dxdxflux + dydyflux})
- contructed_rhs = dtS[subdomain][phase] - div_flux[subdomain][phase]
+ div_flux[subdomain].update({phase: div_flux_x + div_flux_y})
+ contructed_rhs = dtS[subdomain][phase] + div_flux[subdomain][phase]
output['source'][subdomain].update(
{phase: sym.printing.ccode(contructed_rhs)}
)
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
index 13d24f67a13880cf08c2c74295fe8cbaa3661830..511a99b1398426bb8b7384a82ae7a6c28a1c9fc8 100755
--- 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
@@ -19,7 +19,7 @@ datestr = date.strftime("%Y-%m-%d")
# init sympy session
sym.init_printing()
-use_case = "TP-R-2-patch-realistic-pure-dd-new-gravity-rhs"
+use_case = "TP-R-2-patch-realistic-pure-dd"
# solver_tol = 6E-7
max_iter_num = 1000
FEM_Lagrange_degree = 1
@@ -30,26 +30,26 @@ resolutions = {
# 4: 7e-7,
# 8: 7e-7,
# 16: 7e-7,
- 32: 5e-6,
+ 32: 1e-6,
# 64: 7e-7,
# 128: 7e-7,
# 256: 7e-7
}
############ GRID #######################
-timestep_size = 0.001
-number_of_timesteps = 15
+timestep_size = 0.00001
+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
+number_of_timesteps_to_analyse = 4
starttimes = [0.0]
-Lw = 0.04 #/timestep_size
-Lnw= 0.8
+Lw = 0.5 #/timestep_size
+Lnw= 0.5
-lambda_w = 2
-lambda_nw = 8
+lambda_w = 1
+lambda_nw = 1
include_gravity = True
debugflag = True
@@ -72,7 +72,7 @@ else:
write_to_file = {
'space_errornorms': True,
'meshes_and_markers': True,
- 'L_iterations_per_timestep': False,
+ 'L_iterations_per_timestep': True,
'solutions': True,
'absolute_differences': True,
'condition_numbers': analyse_condition,
@@ -190,10 +190,10 @@ lambda_param = {#
}
intrinsic_permeability = {
- 1: {"wetting": 10e-6,
- "nonwetting": 10e-6},
- 2: {"wetting": 10e-6,
- "nonwetting": 10e-6},
+ 1: {"wetting": 1,
+ "nonwetting": 1},
+ 2: {"wetting": 1,
+ "nonwetting": 1},
}
## relative permeabilty functions on subdomain 1
diff --git a/Two-phase-Two-phase/one-patch/TP-one-patch/TP-one-patch-purely-postive-pc.py b/Two-phase-Two-phase/one-patch/TP-one-patch/TP-one-patch-purely-postive-pc.py
index 330fd63e090236ba477ebc792334622ad5ebe271..d0c351558773d9d2a72fed2f0554b964ce38b785 100755
--- a/Two-phase-Two-phase/one-patch/TP-one-patch/TP-one-patch-purely-postive-pc.py
+++ b/Two-phase-Two-phase/one-patch/TP-one-patch/TP-one-patch-purely-postive-pc.py
@@ -39,21 +39,21 @@ resolutions = {
############ GRID #######################
# mesh_resolution = 20
-timestep_size = 0.0005
-number_of_timesteps = 2500
+timestep_size = 0.001
+number_of_timesteps = 25
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 = 0
starttime = 0.0
-Lw = 0.25 #/timestep_size
+Lw = 0.025 #/timestep_size
Lnw=Lw
-lambda_w = 40
-lambda_nw = 40
+lambda_w = 4
+lambda_nw = 4
-include_gravity = False
+include_gravity = True
debugflag = True
analyse_condition = False
diff --git a/Two-phase-Two-phase/two-patch/mesh_studies/TP-TP-2-patch-mesh-study.py b/Two-phase-Two-phase/two-patch/mesh_studies/TP-TP-2-patch-mesh-study.py
index cd8e07900c296d68d2a8284637562d0873142749..bd6ee937074f6c998921da8e032fea90f675c629 100755
--- a/Two-phase-Two-phase/two-patch/mesh_studies/TP-TP-2-patch-mesh-study.py
+++ b/Two-phase-Two-phase/two-patch/mesh_studies/TP-TP-2-patch-mesh-study.py
@@ -23,23 +23,23 @@ use_case = "TP-TP-2-patch-realistic"
# solver_tol = 6E-7
max_iter_num = 300
FEM_Lagrange_degree = 1
-mesh_study = True
+mesh_study = False
resolutions = {
- 1: 1e-6,
- 2: 1e-6,
- 4: 1e-6,
- 8: 1e-6,
- 16: 1e-6,
+ # 1: 1e-6,
+ # 2: 1e-6,
+ # 4: 1e-6,
+ # 8: 1e-6,
+ # 16: 1e-6,
32: 1e-6,
- 64: 1e-6,
- 128: 1e-6,
- 256: 1e-6,
+ # 64: 1e-6,
+ # 128: 1e-6,
+ # 256: 1e-6,
}
############ GRID #######################
# mesh_resolution = 20
timestep_size = 0.001
-number_of_timesteps = 800
+number_of_timesteps = 20
plot_timestep_every = 2
# decide how many timesteps you want analysed. Analysed means, that we write out
# subsequent errors of the L-iteration within the timestep.
@@ -52,8 +52,8 @@ Lnw=Lw
lambda_w = 4
lambda_nw = 4
-include_gravity = False
-debugflag = False
+include_gravity = True
+debugflag = True
analyse_condition = True
if mesh_study: