diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index d5f6413a4b28286fc099a98ec09db244cb48de11..a9d2afc4c958b4e9c1eddc3c42d65a57035d52ea 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -81,7 +81,7 @@ class LDDsimulation(object):
## Private variables
# maximal number of L-iterations that the LDD solver uses.
- self._max_iter_num = 5
+ self._max_iter_num = 100
# TODO rewrite this with regard to the mesh sizes
self.calc_tol = self.tol
# The number of subdomains are counted by self.init_meshes_and_markers()
@@ -180,7 +180,7 @@ class LDDsimulation(object):
self._init_subdomains()
self._init_initial_values()
- def LDDsolver(self, time: float):
+ def LDDsolver(self, time: float, debug: bool = False):
""" calulate the solution on all patches for the next timestep
The heart of the simulation. The function LDDsolver implements the LDD
@@ -215,7 +215,7 @@ class LDDsimulation(object):
subdomain.null_all_interface_iteration_numbers()
# update the time of the source terms before starting the iteration.
# The sources and sinks are the same throughout the iteration.
- self._eval_sources(interpolation_degree = 2,
+ self._eval_sources(interpolation_degree = 3,
time = time,
subdomain_index = ind)
self.update_DirichletBC_dictionary(time = time, subdomain_index = ind)
@@ -228,8 +228,8 @@ class LDDsimulation(object):
subdomain.iteration_number = 0
for phase in subdomain.has_phases:
# set the solution of the previous timestep as new prev_timestep pressure
- subdomain.pressure_prev_timestep[phase].vector().set_local(
- subdomain.pressure[phase].vector().get_local()
+ subdomain.pressure_prev_timestep[phase].assign(
+ subdomain.pressure[phase]
)
### actual iteration starts here
@@ -246,25 +246,32 @@ class LDDsimulation(object):
if not subdomain_calc_isdone[sd_index]:
# set subdomain iteration to current iteration
subdomain.iteration_number = iteration
+ if debug:
+ print(f"Entering iteration {iteration}")
# solve the problem on subdomain
self.Lsolver_step(subdomain_index = sd_index,#
- debug = True,
- )
- # subsequent_iterations_err = self.calc_iteration_error(
- # subdomain_index = sd_index,#
- # error_type = "L2"
- # )
- # # check if error criterion has been met
- # if subsequent_iterations_err < self.calc_tol:
- # subdomain_calc_isdone[sd_index] = True
+ debug = debug,
+ )
+ subsequent_iter_error = dict()
+ for phase in subdomain.has_phases:
+ pa = subdomain.pressure[phase]
+ pa_prev_i = subdomain.pressure_prev_iter[phase]
+ subsequent_iter_error.update(
+ {phase: df.errornorm(pa, pa_prev_i, 'L2')}
+ )
+
+ total_subsequent_iter_error = max(subsequent_iter_error.values())
+ # check if error criterion has been met
+ if total_subsequent_iter_error < subdomain.mesh.hmin()/60:
+ subdomain_calc_isdone[sd_index] = True
# prepare next iteration
# write the newly calculated solution to the inteface dictionaries
# for communication
subdomain.write_pressure_to_interfaces()
for phase in subdomain.has_phases:
- subdomain.pressure_prev_iter[phase].vector().set_local(
- subdomain.pressure[phase].vector().get_local()
+ subdomain.pressure_prev_iter[phase].assign(
+ subdomain.pressure[phase]
)
else:
# one subdomain is done. Check if all subdomains are done to
@@ -356,6 +363,8 @@ class LDDsimulation(object):
if debug and subdomain.mesh.num_cells() < 36:
print("\npressure after solver:\n", subdomain.pressure[phase].vector().get_local())
+
+
## Private methods
def _init_meshes_and_markers(self, subdomain_def_points: tp.List[tp.List[df.Point]] = None,#
mesh_resolution: float = None, #
@@ -609,7 +618,7 @@ class LDDsimulation(object):
pDCa = df.Expression(pDC[phase], domain = mesh, degree = interpolation_degree, t=time)
self.dirichletBC_dfExpression[num].update({phase: pDCa})
else:
- pDCw = self.dirichletBC_dfExpression[num][phase].t = time
+ pDCa = self.dirichletBC_dfExpression[num][phase].t = time
self.outerBC[num].update(
{phase: df.DirichletBC(V[phase], pDCa, boundary_marker, 1)}
@@ -636,9 +645,10 @@ class LDDsimulation(object):
f_expression = {
phase: df.Expression(f, domain = mesh, degree = interpolation_degree, t = time)
}
- # print(f_expression.__dir__())
- # subdomain.source = {'wetting' : df.interpolate(fw, V['wetting'])}
- subdomain.source = f_expression.copy()
+
+ subdomain.source.update(
+ {phase: f_expression[phase]}#
+ )
else:
# note that the default interpolation degree is 2
subdomain.source[phase].t = time
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 03678bd0ccea5ed98ec41702cd42d8ca8148d452..4a4624ead2a738b6dfd30aa99b5bcd714e04a714 100755
--- a/RR-2-patch-test-case/RR-2-patch-test.py
+++ b/RR-2-patch-test-case/RR-2-patch-test.py
@@ -85,7 +85,7 @@ isRichards = {
}
Tmax = 2
-dt1 = 0.01
+dt1 = 0.001
dt2 = dt1
# the timestep is taken as a dictionary to be able to run different timesteps on
# different subdomains.
@@ -185,7 +185,7 @@ write_to_file = {
'L_iterations': True
}
-mesh_resolution = 4
+mesh_resolution = 100
# initialise LDD simulation class
simulation = ldd.LDDsimulation(tol = 1E-12)
@@ -239,7 +239,7 @@ df.File('./subdomain2_interface_marker.pvd') << simulation.subdomain[2].interfac
# df.File('./test_subdomain1_boundary_marker.pvd') << boundary_marker1
# df.File('./test_subdomain2_boundary_marker.pvd') << boundary_marker2
-simulation.LDDsolver(time = 0)
+simulation.LDDsolver(time = 0, debug = True)
# df.info(parameters, True)
# for iterations in range(1):