diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index 85ec4c138c0d94a851d0124fbc7867aed1f18aa3..ede0b8af4d58b4edd256882b06538b15f7ab7b0b 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -100,9 +100,9 @@ class LDDsimulation(object):
## Private variables
# maximal number of L-iterations that the LDD solver uses.
- self._max_iter_num = 5000
+ self._max_iter_num = 1
# TODO rewrite this with regard to the mesh sizes
- self.calc_tol = self.tol
+ # self.calc_tol = self.tol
# list of timesteps that get analyed. Gets initiated by self._init_analyse_timesteps
self.analyse_timesteps = None
# The number of subdomains are counted by self.init_meshes_and_markers()
@@ -143,14 +143,14 @@ class LDDsimulation(object):
########## SOLVER DEFINITION AND PARAMETERS ########
# print("\nLinear Algebra Backends:")
# df.list_linear_algebra_backends()
- # print("\nLinear Solver Methods:")
+ # # print("\nLinear Solver Methods:")
# df.list_linear_solver_methods()
- # print("\nPeconditioners for Krylov Solvers:")
+ # # print("\nPeconditioners for Krylov Solvers:")
# df.list_krylov_solver_preconditioners()
### Define the linear solver to be used.
self.solver = 'bicgstab' #'gmres'#'bicgstab' # biconjugate gradient stabilized method
- self.preconditioner = 'hypre_amg' #'ilu'#'hypre_amg' # incomplete LU factorization
+ self.preconditioner = 'jacobi'#'hypre_amg' #'ilu'#'hypre_amg' # incomplete LU factorization
# dictionary of solver parametrs. This is passed to self._init_subdomains,
# where for each subdomain a sovler object of type self.solver is created
# with these parameters.
@@ -161,6 +161,17 @@ class LDDsimulation(object):
'maximum_iterations': 1000,
'report': False
}
+ # Dictionaries needed by the LDD Solver
+ # dictionary holding bool variables for each subdomain indicating wether
+ # minimal error has been achieved, i.e. the calculation can be considered
+ # done or not.
+ self._calc_isdone_for_subdomain = dict()
+ # dictionary holding for each subdomain a dictionary with bool variables
+ # for each phase indicating wether minimal error has been achieved, i.e.
+ # the calculation can be considered done or not. This dictionary is necessary,
+ # because it can happen, that one phase needs severly less iterations than
+ # the other to achieve the error criterion.
+ self._calc_isdone_for_phase = dict()
def set_parameters(self,#
@@ -285,9 +296,9 @@ class LDDsimulation(object):
# check if the solver should beanalised or not.
if analyse_solver_at_timesteps is not None and np.isin(self.timestep_num, analyse_solver_at_timesteps, assume_unique=True):
print("analyising timestep ...")
- self.LDDsolver(debug = False, analyse_timestep = True)
+ self.LDDsolver(debug = True, analyse_timestep = True)
else:
- self.LDDsolver(debug = False, analyse_timestep = False)
+ self.LDDsolver(debug = True, analyse_timestep = False)
self.t += self.timestep_size
self.timestep_num += 1
@@ -347,18 +358,16 @@ class LDDsimulation(object):
if time is None:
time = self.t
+ error_criterion = 1E-8#min(1E-9, subdomain.mesh.hmax()**2*subdomain.timestep_size)
+
if analyse_timestep:
# dictionary holding filename strings for each subdomain to store
# iterations in, used to analyse the behaviour of the iterates within
# the current time step.
solution_over_iteration_within_timestep = dict()
subsequent_errors_over_iteration = dict()
- # dictionary holding bool variables for each subdomain indicating wether
- # minimal error has been achieved, i.e. the calculation can be considered
- # done or not.
- subdomain_calc_isdone = dict()
# before the iteration gets started all iteration numbers must be nulled
- # and the subdomain_calc_isdone dictionary must be set to False
+ # and the self._calc_isdone_for_subdomain dictionary must be set to False
for ind, subdomain in self.subdomain.items():
# create xdmf files for writing out iterations if analyse_timestep is
# given.
@@ -375,7 +384,7 @@ class LDDsimulation(object):
# subsequent_errors_over_iteration[ind].update( {phase: dict()} )
#
- subdomain_calc_isdone.update({ind: False})
+ self._calc_isdone_for_subdomain.update({ind: False})
# reset all interface[has_interface].current_iteration[ind] = 0, for
# index has_interface in subdomain.has_interface
subdomain.null_all_interface_iteration_numbers()
@@ -386,18 +395,16 @@ class LDDsimulation(object):
self.update_DirichletBC_dictionary(time = time, #
subdomain_index = ind,
)
- # the following only needs to be done when time is not equal 0 since
- # our initialisation functions already took care of setting what follows
- # for the initial iteration step at time = 0.
- if time > self.calc_tol:
- # this is the beginning of the solver, so iteration must be set
- # to 0.
- 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].assign(
- subdomain.pressure[phase]
- )
+ # # this is the beginning of the solver, so iteration must be set
+ # # to 0.
+ 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].assign(
+ subdomain.pressure[phase]
+ )
+ # All isdone flags need to be zero before we start iterating
+ self._calc_isdone_for_phase[ind].update({phase: False})
### actual iteration starts here
# gobal stopping criterion for the iteration.
@@ -410,7 +417,7 @@ class LDDsimulation(object):
for sd_index, subdomain in self.subdomain.items():
# check if the calculation on subdomain has already be marked as
# finished.
- if not subdomain_calc_isdone[sd_index]:
+ if not self._calc_isdone_for_subdomain[sd_index]:
# set subdomain iteration to current iteration
subdomain.iteration_number = iteration
if debug:
@@ -430,6 +437,10 @@ class LDDsimulation(object):
subsequent_iter_error.update(
{phase: error_calculated}
)
+ # set flag to stop solving for that phase if error_criterion is met.
+ if subsequent_iter_error[phase] < error_criterion:
+ self._calc_isdone_for_phase[sd_index].update({phase: True})
+
if debug:
print(f"the subsequent error in iteration {iteration} for phase {phase} is {subsequent_iter_error[phase]}")
@@ -458,16 +469,18 @@ class LDDsimulation(object):
print("max(subsequent_iter_error.values()): ", max(subsequent_iter_error.values()))
total_subsequent_iter_error = max(subsequent_iter_error.values())
# check if error criterion has been met
- error_criterion = 1E-9#min(1E-9, subdomain.mesh.hmax()**2*subdomain.timestep_size)
+
if debug:
print(f" total_subsequent_error in iter {iteration} = {total_subsequent_iter_error }\n")
print(f"the maximum distance between any to points in a cell mesh.hmax() on subdomain{sd_index} is h={subdomain.mesh.hmax()}")
print(f"the minimal distance between any to points in a cell mesh.hmin() on subdomain{sd_index} is h={subdomain.mesh.hmin()}")
print(f"the error criterion is tol = {error_criterion}")
- if total_subsequent_iter_error < error_criterion:
+ wetting_done = self._calc_isdone_for_phase[sd_index]['wetting']
+ nonwetting_done = self._calc_isdone_for_phase[sd_index]['nonwetting']
+ if nonwetting_done and wetting_done:
if debug:
print(f"subdomain{sd_index} reachd error tol = {error_criterion} in iteration {iteration}.")
- subdomain_calc_isdone[sd_index] = True
+ self._calc_isdone_for_subdomain[sd_index] = True
# prepare next iteration
# write the newly calculated solution to the inteface dictionaries
@@ -491,7 +504,7 @@ class LDDsimulation(object):
all_subdomains_are_done = True
# then check if all domains are done. If not, reset
# all_subdomains_are_done to False.
- for subdomain, isdone in subdomain_calc_isdone.items():
+ for subdomain, isdone in self._calc_isdone_for_subdomain.items():
if not isdone:
all_subdomains_are_done = False
#TODO check if maybe I still need to do
@@ -522,61 +535,65 @@ class LDDsimulation(object):
to write out each iteration.
"""
subdomain = self.subdomain[subdomain_index]
+ calc_isdone_for_phase = self._calc_isdone_for_phase[subdomain_index]
iteration = subdomain.iteration_number
# calculate the gli terms on the right hand side and store
subdomain.calc_gli_term()
for phase in subdomain.has_phases:
- # extract L-scheme form and rhs (without gli term) from subdomain.
- governing_problem = subdomain.governing_problem(phase = phase)
- form = governing_problem['form']
- rhs_without_gli = governing_problem['rhs_without_gli']
- # assemble the form and rhs
- form_assembled = df.assemble(form)
- rhs_without_gli_assembled = df.assemble(rhs_without_gli)
- # access the assembled gli term for the rhs.
- gli_term_assembled = subdomain.gli_assembled[phase]
- # subdomain.calc_gli_term() asslembles gli but on the left hand side
- # so gli_term_assembled needs to actually be subtracted from the rhs.
- if debug and subdomain.mesh.num_cells() < 36:
- print("\nSystem before applying outer boundary conditions:")
- print(f"phase = {phase}: rhs_without_gli:\n", rhs_without_gli_assembled.get_local())
- rhs_without_gli_assembled.set_local(rhs_without_gli_assembled.get_local() - gli_term_assembled)
- rhs_assembled = rhs_without_gli_assembled
-
- if debug and subdomain.mesh.num_cells() < 36:
- # print(f"phase = {phase}: form_assembled:\n", form_assembled.array())
- # print(f"phase = {phase}: rhs_without_gli:\n", rhs_without_gli_assembled.get_local())
- print(f"phase = {phase}: gli_term:\n", gli_term_assembled)
- print(f"phase = {phase}: rhs_assembled:\n", rhs_assembled.get_local())
-
- # apply outer Dirichlet boundary conditions if present.
- if subdomain.outer_boundary is not None:
- # apply boundary dirichlet conditions for all outer boundaries.
- for index, boundary in enumerate(subdomain.outer_boundary):
- self.outerBC[subdomain_index][index][phase].apply(form_assembled, rhs_assembled)
-
- if debug and subdomain.mesh.num_cells() < 36:
- print("\nSystem after applying outer boundary conditions:")
- # print(f"phase = {phase}: form_assembled:\n", form_assembled.array())
- print(f"phase = {phase}: rhs_assembled:\n", rhs_assembled.get_local())
-
- # # set previous iteration as initial guess for the linear solver
- # pi_prev = subdomain.pressure_prev_iter[phase].vector().get_local()
- # subdomain.pressure[phase].vector().set_local(pi_prev)
-
- if debug and subdomain.mesh.num_cells() < 36:
- print("\npressure before solver:\n", subdomain.pressure[phase].vector().get_local())
-
- linear_solver = subdomain.linear_solver
- if linear_solver is None:
- df.solve(form_assembled, subdomain.pressure[phase].vector(), rhs_assembled, self.solver, self.preconditioner)
- else:
- linear_solver.solve(form_assembled, subdomain.pressure[phase].vector(), rhs_assembled)
-
- if debug and subdomain.mesh.num_cells() < 36:
- print("\npressure after solver:\n", subdomain.pressure[phase].vector().get_local())
+ if not calc_isdone_for_phase[phase]:
+ # extract L-scheme form and rhs (without gli term) from subdomain.
+ governing_problem = subdomain.governing_problem(phase = phase)
+ form = governing_problem['form']
+ rhs_without_gli = governing_problem['rhs_without_gli']
+ # assemble the form and rhs
+ form_assembled = df.assemble(form)
+ rhs_without_gli_assembled = df.assemble(rhs_without_gli)
+ # access the assembled gli term for the rhs.
+ gli_term_assembled = subdomain.gli_assembled[phase]
+ # subdomain.calc_gli_term() asslembles gli but on the left hand side
+ # so gli_term_assembled needs to actually be subtracted from the rhs.
+ if debug and subdomain.mesh.num_cells() < 36:
+ print("\nSystem before applying outer boundary conditions:")
+ print(f"phase = {phase}: rhs_without_gli:\n", rhs_without_gli_assembled.get_local())
+ rhs_without_gli_assembled.set_local(rhs_without_gli_assembled.get_local() - gli_term_assembled)
+ rhs_assembled = rhs_without_gli_assembled
+
+ if debug and subdomain.mesh.num_cells() < 36:
+ # print(f"phase = {phase}: form_assembled:\n", form_assembled.array())
+ # print(f"phase = {phase}: rhs_without_gli:\n", rhs_without_gli_assembled.get_local())
+ print(f"phase = {phase}: gli_term:\n", gli_term_assembled)
+ print(f"phase = {phase}: rhs_assembled:\n", rhs_assembled.get_local())
+
+ # apply outer Dirichlet boundary conditions if present.
+ if subdomain.outer_boundary is not None:
+ # apply boundary dirichlet conditions for all outer boundaries.
+ for index, boundary in enumerate(subdomain.outer_boundary):
+ self.outerBC[subdomain_index][index][phase].apply(form_assembled, rhs_assembled)
+
+ if debug and subdomain.mesh.num_cells() < 36:
+ print("\nSystem after applying outer boundary conditions:")
+ # print(f"phase = {phase}: form_assembled:\n", form_assembled.array())
+ print(f"phase = {phase}: rhs_assembled:\n", rhs_assembled.get_local())
+
+ # # set previous iteration as initial guess for the linear solver
+ # pi_prev = subdomain.pressure_prev_iter[phase].vector().get_local()
+ # subdomain.pressure[phase].vector().set_local(pi_prev)
+
+ if debug and subdomain.mesh.num_cells() < 36:
+ print("\npressure before solver:\n", subdomain.pressure[phase].vector().get_local())
+
+ linear_solver = subdomain.linear_solver
+ if linear_solver is None:
+ df.solve(form_assembled, subdomain.pressure[phase].vector(), rhs_assembled, self.solver, self.preconditioner)
+ else:
+ linear_solver.solve(form_assembled, subdomain.pressure[phase].vector(), rhs_assembled)
+ if debug and subdomain.mesh.num_cells() < 36:
+ print("\npressure after solver:\n", subdomain.pressure[phase].vector().get_local())
+ else:
+ print(f"calculation for phase {phase} is already done in timestep = {self.timestep_num}, t = {self.t}. skipping ..." )
+ pass
def _init_solution_files(self):
""" set up solution files for saving the solution of the LDD scheme for
@@ -588,7 +605,7 @@ class LDDsimulation(object):
for subdom_ind, subdomain in self.subdomain.items():
tau = subdomain.timestep_size
L = subdomain.L
- h = subdomain.mesh.hmin()
+ h = subdomain.mesh.hmax()
if subdomain.isRichards:
Lam = subdomain.lambda_param['wetting']
name1 ="subdomain" + "{}".format(subdom_ind)
@@ -836,7 +853,7 @@ class LDDsimulation(object):
# marker value gets internally set to global_index+1.
# The reason is that we want to mark outer boundaries with the value 0.
self.interface[num].mark(self.interface_marker, self.interface[num].marker_value)
- print("Global interface vertices:")
+ # print("Global interface vertices:")
self.interface[num].init_interface_values(self.interface_marker, self.interface[num].marker_value)
if self.write2file['meshes_and_markers']:
@@ -887,6 +904,8 @@ class LDDsimulation(object):
# for parameter, value in self.linear_solver.parameters.items():
# print(f"parameter: {parameter} = {value}")
# df.info(solver_param, True)
+ # setup self._calc_isdone_for_phase dictionary.
+ self._calc_isdone_for_phase.update({subdom_num: dict()})
@@ -1059,21 +1078,15 @@ class LDDsimulation(object):
subdomain = self.subdomain[subdomain_index]
for phase in subdomain.has_phases:
- if np.abs(time-self.starttime) < self.calc_tol:
+ if np.abs(time-self.starttime) < self.tol:
# here t = 0 and we have to initialise the sources.
- mesh = subdomain.mesh
- # p0 contains both pw_0 and pnw_0 for subdomain[subdomain_index]
f = self.sources[subdomain_index][phase]
- # note that the default interpolation degree is 2
- f_expression = {
- phase: df.Expression(f, domain = mesh, degree = interpolation_degree, t = time)
- }
-
subdomain.source.update(
- {phase: f_expression[phase]}#
+ {phase: df.Expression(f, domain = subdomain.mesh,#
+ degree = interpolation_degree, #
+ t = time)}
)
else:
- # note that the default interpolation degree is 2
subdomain.source[phase].t = time
def _add_parameters_to_output_dirname(self):