diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index 0b5678d30b85b631fd5bb2cdbf1e9313b1929fa1..a9ddac2ccfe46c2c93851f5f8050a5ab11c27df2 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -162,7 +162,7 @@ class LDDsimulation(object):
'absolute_tolerance': 1E-12,
'relative_tolerance': 1E-10,
'maximum_iterations': 1000,
- 'report': True
+ 'report': False
}
self.debug_solver = debug
@@ -897,6 +897,7 @@ class LDDsimulation(object):
self.interface[num].mark(self.interface_marker, self.interface[num].marker_value)
# print("Global interface vertices:")
self.interface[num].init_interface_values(self.interface_marker, self.interface[num].marker_value)
+ # print(self.interface[num])
if self.write2file['meshes_and_markers']:
filepath = self.output_dir+"interface_marker_global.pvd"
@@ -954,6 +955,7 @@ class LDDsimulation(object):
# for parameter, value in self.linear_solver.parameters.items():
# print(f"parameter: {parameter} = {value}")
# df.info(solver_param, True)
+ print(self.subdomain[subdom_num])
def _get_interfaces(self, subdomain_index: int, #
diff --git a/LDDsimulation/boundary_and_interface.py b/LDDsimulation/boundary_and_interface.py
index 55250d3bc922ce89f8788e3e42078e1fbf736411..12b3ec2b276d1252bfa8d12f7dc1103c9a6e6eb7 100644
--- a/LDDsimulation/boundary_and_interface.py
+++ b/LDDsimulation/boundary_and_interface.py
@@ -1,3 +1,4 @@
+#!/bin/python3
""" BoundaryPart and Interface class
This module defines the classes BoundaryPart and its child, the class Interface
@@ -262,6 +263,15 @@ class Interface(BoundaryPart):
# make the parent class methods available.
super().__init__(**kwds)
+ # MAGIC METHODS
+ def __str__(self):
+ """print out a few usefull informations of the interface"""
+ string = f"\nInterface{self.global_index} with "\
+ +f" global marker_value = {self.marker_value}.\n"\
+ +f" Neighbouring subdomains are subdomain{self.adjacent_subdomains[0]}"\
+ +f" and subdomain{self.adjacent_subdomains[1]}"
+ return string
+
#### Public Methods #####
# def set_domain_index(self, domain_index: int):
# self.domain_index = domain_index
diff --git a/LDDsimulation/domainPatch.py b/LDDsimulation/domainPatch.py
index c0aad7b3c9d5ded1981304e9ae6312e6c2289914..94fa2ca690dc12e7adf69f2614886643bf8ef8b8 100644
--- a/LDDsimulation/domainPatch.py
+++ b/LDDsimulation/domainPatch.py
@@ -234,6 +234,30 @@ class DomainPatch(df.SubDomain):
# END constructor
+
+ # MAGIC METHODS
+ def __str__(self):
+ """print out a few usefull informations of the subdomain"""
+ model = "two-phase"
+ if self.isRichards:
+ model = "Richards"
+
+ interface_string = ["\nMy interfaces are:"]
+ for glob_if_index in self.has_interface:
+ interface = self.interface[glob_if_index]
+ interface_string.append(interface.__str__())
+ common_int_dof_str = ".\n Dofs common with other interfaces:"\
+ +"{}".format([(ind, dof) for ind, dof in self._interface_has_common_dof_indices[glob_if_index].items()])
+ interface_string.append(common_int_dof_str)
+
+ interface_string = " ".join(interface_string)
+ string = f"\nI'm subdomain{self.subdomain_index} and assume {model} model."\
+ +f"\nI have the interfaces: {self.has_interface} and "\
+ +f"outer boundaries: {[ind for ind in self.outer_boundary_def_points.keys()]}.\n"\
+ + interface_string
+ return string
+
+
#### PUBLIC METHODS
def write_pressure_to_interfaces(self, previous_iter: int = False):
""" save the interface values of self.pressure to all neighbouring interfaces
@@ -422,7 +446,8 @@ class DomainPatch(df.SubDomain):
pa_prev = self.pressure_prev_timestep[phase]
# testfunction
phi_a = self.testfunction[phase]
- z_a = self.gravity_term[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
@@ -482,8 +507,11 @@ class DomainPatch(df.SubDomain):
# dictionaries.
for phase in self.has_phases:
for interf_ind in self.has_interface:
- # self._calc_gli_term_on(interface, iteration)
interface = self.interface[interf_ind]
+ neighbour = interface.neighbour[subdomain]
+ neighbour_iter_num = interface.current_iteration[neighbour]
+ # needed for the read_gli_dofs() functions
+ interface_dofs = self._dof_indices_of_interface[interf_ind]
if debug:
print(f"Interface{interf_ind}.gli_term_prev[{subdomain}][{phase}]=\n",interface.gli_term_prev[subdomain][phase])
print(f"before writing to interface[{interf_ind}] the values gli_assembled[{phase}]=\n",gli_assembled_tmp[phase])
@@ -630,6 +658,11 @@ class DomainPatch(df.SubDomain):
for ind in self.has_interface:
# self._calc_gli_term_on(interface, iteration)
interface = self.interface[ind]
+ neighbour = interface.neighbour[subdomain]
+ neighbour_iter_num = interface.current_iteration[neighbour]
+ # needed for the read_gli_dofs() functions
+ interface_dofs = self._dof_indices_of_interface[ind]
+
if debug:
print(f"Interface{ind}.gli_term[{subdomain}][{phase}]=\n",interface.gli_term[subdomain][phase])
print(f"before writing to interface[{ind}] gli_assembled[{phase}]=\n",self.gli_assembled[phase])
@@ -669,6 +702,8 @@ class DomainPatch(df.SubDomain):
# if we don't save it to gli_term_prev of the neighbour.
# similarly for the pressure values.
if phase == 'wetting' or not interface.neighbour_isRichards[subdomain]:
+ # print(f"we are on subdomain{subdomain} and interface{interface.global_index} (alt index:{ind}) has neighbouring \n",
+ # f"subdomains: subdomain{subdomain} and subdomain{neighbour}")
interface.gli_term_prev[neighbour].update(#
{phase: interface.gli_term[neighbour][phase].copy()}
)
diff --git a/RR-2-patch-test-case/RR-2-patch-test.py b/RR-2-patch-test-case/RR-2-patch-test.py
index 4cc538c278734e90d43c1c4fb14b96e0ca514f10..f3d946e7b7de6046c95548681e567ca5727eea92 100755
--- a/RR-2-patch-test-case/RR-2-patch-test.py
+++ b/RR-2-patch-test-case/RR-2-patch-test.py
@@ -11,13 +11,13 @@ import functools as ft
##### Domain and Interface ####
# global simulation domain domain
-sub_domain0_vertices = [df.Point(0.0,0.0), #
- df.Point(1.0,0.0),#
+sub_domain0_vertices = [df.Point(0.0,-1.0), #
+ df.Point(1.0,-1.0),#
df.Point(1.0,1.0),#
df.Point(0.0,1.0)]
# interface between subdomain1 and subdomain2
-interface12_vertices = [df.Point(0.0, 0.5),
- df.Point(1.0, 0.5) ]
+interface12_vertices = [df.Point(0.0, 0.0),
+ df.Point(1.0, 0.0) ]
# subdomain1.
sub_domain1_vertices = [interface12_vertices[0],
interface12_vertices[1],
@@ -29,22 +29,22 @@ sub_domain1_vertices = [interface12_vertices[0],
# the Dirichlet boundary conditions. If a domain is completely internal, the
# dictionary entry should be 0: None
subdomain1_outer_boundary_verts = {
- 0: [interface12_vertices[0], #
- df.Point(0.0,1.0), #
+ 0: [interface12_vertices[1], #
df.Point(1.0,1.0), #
- interface12_vertices[1]]
+ df.Point(0.0,1.0), #
+ interface12_vertices[0]]
}
# subdomain2
-sub_domain2_vertices = [df.Point(0.0,0.0),
- df.Point(1.0,0.0),
+sub_domain2_vertices = [df.Point(0.0,-1.0),
+ df.Point(1.0,-1.0),
interface12_vertices[1],
interface12_vertices[0] ]
subdomain2_outer_boundary_verts = {
- 0: [interface12_vertices[1], #
- df.Point(1.0,0.0), #
- df.Point(0.0,0.0), #
- interface12_vertices[0]]
+ 0: [interface12_vertices[0], #
+ df.Point(0.0,-1.0), #
+ df.Point(1.0,-1.0), #
+ interface12_vertices[1]]
}
# subdomain2_outer_boundary_verts = {
# 0: [interface12_vertices[0], df.Point(0.0,0.0)],#
@@ -85,10 +85,10 @@ isRichards = {
}
##################### MESH #####################################################
-mesh_resolution = 100
+mesh_resolution = 5
##################### TIME #####################################################
-timestep_size = 4*0.0001
-number_of_timesteps = 3000
+timestep_size = 0.01
+number_of_timesteps = 300
# 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
@@ -101,6 +101,13 @@ viscosity = {#
2 : {'wetting' :1}
}
+densities = {#
+# subdom_num : viscosity
+ 1 : {'wetting' :1}, #
+ 2 : {'wetting' :1}
+}
+
+
porosity = {#
# subdom_num : porosity
1 : 1,#
@@ -115,8 +122,8 @@ L = {#
lambda_param = {#
# subdom_num : lambda parameter for the L-scheme
- 1 : {'wetting' :100},#
- 2 : {'wetting' :100}
+ 1 : {'wetting' :40},#
+ 2 : {'wetting' :40}
}
## relative permeabilty functions on subdomain 1
@@ -163,17 +170,17 @@ sat_pressure_relationship = {#
source_expression = {
1: {'wetting': '4.0/pow(1 + x[0]*x[0] + x[1]*x[1], 2) - t/sqrt( pow(1 + t*t, 3)*(1 + x[0]*x[0] + x[1]*x[1]) )'},
- 2: {'wetting': '2.0*(1-x[1]*x[1])/pow(1 + x[1]*x[1], 2) - 2*t/(3*pow( pow(1 + t*t, 4)*(1 + x[1]*x[1]), 1/3))'}
+ 2: {'wetting': '2.0*(1-x[0]*x[0])/pow(1 + x[0]*x[0], 2) - 2*t/(3*pow( pow(1 + t*t, 4)*(1 + x[0]*x[0]), 1/3))'}
}
initial_condition = {
1: {'wetting': '-(x[0]*x[0] + x[1]*x[1])'},#
- 2: {'wetting': '-x[1]*x[1]'}
+ 2: {'wetting': '-x[0]*x[0]'}
}
exact_solution = {
1: {'wetting': '1.0 - (1.0 + t*t)*(1.0 + x[0]*x[0] + x[1]*x[1])'},#
- 2: {'wetting': '1.0 - (1.0 + t*t)*(1.0 + x[1]*x[1])'}
+ 2: {'wetting': '1.0 - (1.0 + t*t)*(1.0 + x[0]*x[0])'}
}
# similary to the outer boundary dictionary, if a patch has no outer boundary
@@ -186,8 +193,8 @@ exact_solution = {
# the actual expression needed for the dirichlet condition for both phases if present.
dirichletBC = {
#subdomain index: {outer boudary part index: {phase: expression}}
- 1: { 0: {'wetting': '1.0 - (1.0 + t*t)*(1.0 + x[0]*x[0] + x[1]*x[1])'}},
- 2: { 0: {'wetting': '1.0 - (1.0 + t*t)*(1.0 + x[1]*x[1])'}}
+ 1: { 0: {'wetting': exact_solution[1]['wetting']}},
+ 2: { 0: {'wetting': exact_solution[2]['wetting']}}
}
# def saturation(pressure, subdomain_index):
@@ -202,7 +209,7 @@ write_to_file = {
}
# initialise LDD simulation class
-simulation = ldd.LDDsimulation(tol = 1E-14)
+simulation = ldd.LDDsimulation(tol = 1E-14, debug = False, LDDsolver_tol=1E-9)
simulation.set_parameters(output_dir = "./output/",#
subdomain_def_points = subdomain_def_points,#
isRichards = isRichards,#
@@ -224,10 +231,13 @@ simulation.set_parameters(output_dir = "./output/",#
initial_conditions = initial_condition,#
dirichletBC_expression_strings = dirichletBC,#
exact_solution = exact_solution,#
+ densities = densities,#
+ include_gravity = False,#
write2file = write_to_file,#
)
simulation.initialise()
+print(simulation.__dict__)
simulation.run()
# simulation.LDDsolver(time = 0, debug = True, analyse_timestep = True)
# df.info(parameters, True)
diff --git a/TP-TP-patch-test-case/TP-TP-2-patch-test.py b/TP-TP-patch-test-case/TP-TP-2-patch-test.py
index daf7b639020bc9d89b2b42d4d3d9218ad1cbc68f..222ab56bdf1c783c62177dcb786a512a1ba8ccad 100755
--- a/TP-TP-patch-test-case/TP-TP-2-patch-test.py
+++ b/TP-TP-patch-test-case/TP-TP-2-patch-test.py
@@ -88,10 +88,10 @@ isRichards = {
}
-############ GRID ########################ü
+############ GRID #######################ü
mesh_resolution = 20
-timestep_size = 0.0001
-number_of_timesteps = 100
+timestep_size = 0.001
+number_of_timesteps = 3000
# 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
@@ -100,9 +100,9 @@ starttime = 0
viscosity = {#
# subdom_num : viscosity
1 : {'wetting' :1,
- 'nonwetting': 1/50}, #
+ 'nonwetting': 1}, #
2 : {'wetting' :1,
- 'nonwetting': 1/50}
+ 'nonwetting': 1}
}
porosity = {#
@@ -113,18 +113,19 @@ porosity = {#
L = {#
# subdom_num : subdomain L for L-scheme
- 1 : {'wetting' :0.4,
- 'nonwetting': 0.8},#
- 2 : {'wetting' :0.4,
- 'nonwetting': 0.8}
+ 1 : {'wetting' :0.25,
+ 'nonwetting': 0.25},#
+ 2 : {'wetting' :0.25,
+ 'nonwetting': 0.25}
}
+l_param = 40
lambda_param = {#
# subdom_num : lambda parameter for the L-scheme
- 1 : {'wetting' :10,
- 'nonwetting': 100},#
- 2 : {'wetting' :10,
- 'nonwetting': 100}
+ 1 : {'wetting' :l_param,
+ 'nonwetting': l_param},#
+ 2 : {'wetting' :l_param,
+ 'nonwetting': l_param}
}
## relative permeabilty functions on subdomain 1
@@ -134,7 +135,7 @@ def rel_perm1w(s):
def rel_perm1nw(s):
# relative permeabilty nonwetting on subdomain1
- return (1-s)**3
+ return (1-s)**2
_rel_perm1w = ft.partial(rel_perm1w)
_rel_perm1nw = ft.partial(rel_perm1nw)
@@ -145,7 +146,7 @@ subdomain1_rel_perm = {
## relative permeabilty functions on subdomain 2
def rel_perm2w(s):
# relative permeabilty wetting on subdomain2
- return s**3
+ return s**2
def rel_perm2nw(s):
# relative permeabilty nonwetting on subdosym.cos(0.8*t - (0.8*x + 1/7*y))main2
return (1-s)**2
@@ -194,20 +195,20 @@ x, y = sym.symbols('x[0], x[1]') # needed by UFL
t = sym.symbols('t', positive=True)
#f = -sym.diff(u, x, 2) - sym.diff(u, y, 2) # -Laplace(u)
#f = sym.simplify(f) # simplify f
-p1_w = 1 - (1+t)*(1 + x**2 + y**2)
-p1_nw = t*(1-y)**2 - sym.sqrt(2+t**2)*(1-y)
+p1_w = 1 - (1+t**2)*(1 + x**2 + (y-0.5)**2)
+p1_nw = t*(1-(y-0.5) - x**2)**2 - sym.sqrt(2+t**2)*(1-(y-0.5))
#dtS1_w = sym.diff(S_pc_rel_sym[1](p1_nw - p1_w), t, 1)
#dtS1_nw = -sym.diff(S_pc_rel_sym[1](p1_nw - p1_w), t, 1)
-dtS1_w = porosity[1]*sym.diff(S_pc_rel[1](p1_nw - p1_w), t, 1)
-dtS1_nw = -porosity[1]*sym.diff(S_pc_rel[1](p1_nw - p1_w), t, 1)
+dtS1_w = porosity[1]*sym.diff(sym.Piecewise((S_pc_rel[1](p1_nw - p1_w), (p1_nw - p1_w) > 0), (1, True) ), t, 1)
+dtS1_nw = -porosity[1]*sym.diff(sym.Piecewise((S_pc_rel[1](p1_nw - p1_w), (p1_nw - p1_w) > 0), (1, True) ), t, 1)
print("dtS1_w = ", dtS1_w, "\n")
print("dtS1_nw = ", dtS1_nw, "\n")
#dxdxflux1_w = -sym.diff(relative_permeability[1]['wetting'](S_pc_rel_sym[1](p1_nw - p1_w))*sym.diff(p1_w, x, 1), x, 1)
#dydyflux1_w = -sym.diff(relative_permeability[1]['wetting'](S_pc_rel_sym[1](p1_nw - p1_w))*sym.diff(p1_w, y, 1), y, 1)
-dxdxflux1_w = -1/viscosity[1]['wetting']*sym.diff(relative_permeability[1]['wetting'](S_pc_rel[1](p1_nw - p1_w))*sym.diff(p1_w, x, 1), x, 1)
-dydyflux1_w = -1/viscosity[1]['wetting']*sym.diff(relative_permeability[1]['wetting'](S_pc_rel[1](p1_nw - p1_w))*sym.diff(p1_w, y, 1), y, 1)
+dxdxflux1_w = -1/viscosity[1]['wetting']*sym.diff(relative_permeability[1]['wetting'](sym.Piecewise((S_pc_rel[1](p1_nw - p1_w), (p1_nw - p1_w) > 0), (1, True) ))*sym.diff(p1_w, x, 1), x, 1)
+dydyflux1_w = -1/viscosity[1]['wetting']*sym.diff(relative_permeability[1]['wetting'](sym.Piecewise((S_pc_rel[1](p1_nw - p1_w), (p1_nw - p1_w) > 0), (1, True) ))*sym.diff(p1_w, y, 1), y, 1)
rhs1_w = dtS1_w + dxdxflux1_w + dydyflux1_w
rhs1_w = sym.printing.ccode(rhs1_w)
@@ -219,28 +220,28 @@ print("rhs_w = ", rhs1_w, "\n")
#dxdxflux1_nw = -sym.diff(relative_permeability[1]['nonwetting'](S_pc_rel_sym[1](p1_nw - p1_w))*sym.diff(p1_nw, x, 1), x, 1)
#dydyflux1_nw = -sym.diff(relative_permeability[1]['nonwetting'](S_pc_rel_sym[1](p1_nw - p1_w))*sym.diff(p1_nw, y, 1), y, 1)
-dxdxflux1_nw = -1/viscosity[1]['nonwetting']*sym.diff(relative_permeability[1]['nonwetting'](1-S_pc_rel[1](p1_nw - p1_w))*sym.diff(p1_nw, x, 1), x, 1)
-dydyflux1_nw = -1/viscosity[1]['nonwetting']*sym.diff(relative_permeability[1]['nonwetting'](1-S_pc_rel[1](p1_nw - p1_w))*sym.diff(p1_nw, y, 1), y, 1)
+dxdxflux1_nw = -1/viscosity[1]['nonwetting']*sym.diff(relative_permeability[1]['nonwetting'](1-sym.Piecewise((S_pc_rel[1](p1_nw - p1_w), (p1_nw - p1_w) > 0), (1, True) ))*sym.diff(p1_nw, x, 1), x, 1)
+dydyflux1_nw = -1/viscosity[1]['nonwetting']*sym.diff(relative_permeability[1]['nonwetting'](1-sym.Piecewise((S_pc_rel[1](p1_nw - p1_w), (p1_nw - p1_w) > 0), (1, True) ))*sym.diff(p1_nw, y, 1), y, 1)
rhs1_nw = dtS1_nw + dxdxflux1_nw + dydyflux1_nw
rhs1_nw = sym.printing.ccode(rhs1_nw)
print("rhs_nw = ", rhs1_nw, "\n")
## subdomain2
-p2_w = 1 - (1+t)*(1 + x**2 + 0.25*y)
-p2_nw = t*(1-y)**2 - sym.sqrt(2+t**2)*(1-y) +(0.5-y)*(t/(0.1+x))
+p2_w = 1 - (1+t**2)*(1 + x**2)
+p2_nw = t*(1- x**2)**2 - sym.sqrt(2+t**2)*(1-(y-0.5))
#dtS2_w = sym.diff(S_pc_rel_sym[2](p2_nw - p2_w), t, 1)
#dtS2_nw = -sym.diff(S_pc_rel_sym[2](p2_nw - p2_w), t, 1)
-dtS2_w = porosity[2]*sym.diff(S_pc_rel[2](p2_nw - p2_w), t, 1)
-dtS2_nw = -porosity[2]*sym.diff(S_pc_rel[2](p2_nw - p2_w), t, 1)
+dtS2_w = porosity[2]*sym.diff(sym.Piecewise((sym.Piecewise((S_pc_rel[2](p2_nw - p2_w), (p2_nw - p2_w) > 0), (1, True) ), (p2_nw - p2_w) > 0), (1, True) ), t, 1)
+dtS2_nw = -porosity[2]*sym.diff(sym.Piecewise((S_pc_rel[2](p2_nw - p2_w), (p2_nw - p2_w) > 0), (1, True) ), t, 1)
print("dtS2_w = ", dtS2_w, "\n")
print("dtS2_nw = ", dtS2_nw, "\n")
#dxdxflux2_w = -sym.diff(relative_permeability[2]['wetting'](S_pc_rel_sym[2](p2_nw - p2_w))*sym.diff(p2_w, x, 1), x, 1)
#dydyflux2_w = -sym.diff(relative_permeability[2]['wetting'](S_pc_rel_sym[2](p2_nw - p2_w))*sym.diff(p2_w, y, 1), y, 1)
-dxdxflux2_w = -1/viscosity[2]['wetting']*sym.diff(relative_permeability[2]['wetting'](S_pc_rel[2](p2_nw - p2_w))*sym.diff(p2_w, x, 1), x, 1)
-dydyflux2_w = -1/viscosity[2]['wetting']*sym.diff(relative_permeability[2]['wetting'](S_pc_rel[2](p2_nw - p2_w))*sym.diff(p2_w, y, 1), y, 1)
+dxdxflux2_w = -1/viscosity[2]['wetting']*sym.diff(relative_permeability[2]['wetting'](sym.Piecewise((S_pc_rel[2](p2_nw - p2_w), (p2_nw - p2_w) > 0), (1, True) ))*sym.diff(p2_w, x, 1), x, 1)
+dydyflux2_w = -1/viscosity[2]['wetting']*sym.diff(relative_permeability[2]['wetting'](sym.Piecewise((S_pc_rel[2](p2_nw - p2_w), (p2_nw - p2_w) > 0), (1, True) ))*sym.diff(p2_w, y, 1), y, 1)
rhs2_w = dtS2_w + dxdxflux2_w + dydyflux2_w
rhs2_w = sym.printing.ccode(rhs2_w)
@@ -252,8 +253,8 @@ print("rhs2_w = ", rhs2_w, "\n")
#dxdxflux2_nw = -sym.diff(relative_permeability[2]['nonwetting'](S_pc_rel_sym[2](p2_nw - p2_w))*sym.diff(p2_nw, x, 1), x, 1)
#dydyflux2_nw = -sym.diff(relative_permeability[2]['nonwetting'](S_pc_rel_sym[2](p2_nw - p2_w))*sym.diff(p2_nw, y, 1), y, 1)
-dxdxflux2_nw = -1/viscosity[2]['nonwetting']*sym.diff(relative_permeability[2]['nonwetting'](1-S_pc_rel[2](p2_nw - p2_w))*sym.diff(p2_nw, x, 1), x, 1)
-dydyflux2_nw = -1/viscosity[2]['nonwetting']*sym.diff(relative_permeability[2]['nonwetting'](1-S_pc_rel[2](p2_nw - p2_w))*sym.diff(p2_nw, y, 1), y, 1)
+dxdxflux2_nw = -1/viscosity[2]['nonwetting']*sym.diff(relative_permeability[2]['nonwetting'](1-sym.Piecewise((S_pc_rel[2](p2_nw - p2_w), (p2_nw - p2_w) > 0), (1, True) ))*sym.diff(p2_nw, x, 1), x, 1)
+dydyflux2_nw = -1/viscosity[2]['nonwetting']*sym.diff(relative_permeability[2]['nonwetting'](1-sym.Piecewise((S_pc_rel[2](p2_nw - p2_w), (p2_nw - p2_w) > 0), (1, True) ))*sym.diff(p2_nw, y, 1), y, 1)
rhs2_nw = dtS2_nw + dxdxflux2_nw + dydyflux2_nw
rhs2_nw = sym.printing.ccode(rhs2_nw)
@@ -318,7 +319,7 @@ write_to_file = {
# initialise LDD simulation class
-simulation = ldd.LDDsimulation(tol = 1E-14)
+simulation = ldd.LDDsimulation(tol = 1E-14, LDDsolver_tol = 1E-9, debug = False)
simulation.set_parameters(output_dir = "./output/",#
subdomain_def_points = subdomain_def_points,#
isRichards = isRichards,#