diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index c38c71c91b8e959aee88d37afd9c07b9ba3b860d..96681d1d5e12c039c913ad40e97ed2914541b1ba 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -34,7 +34,8 @@ class LDDsimulation(object):
LDDsolver_tol: float = 1E-8,
# maximal number of L-iterations that the LDD solver uses.
max_iter_num: int = 1000,
- FEM_Lagrange_degree: int = 1):
+ FEM_Lagrange_degree: int = 1,
+ plot_timestep_every: int = 1):
hlp.print_once("\n ###############################################\n",
"############# Begin Simulation ################\n",
"###############################################\n")
@@ -115,6 +116,10 @@ class LDDsimulation(object):
# self.calc_tol = self.tol
# list of timesteps that get analyed. Gets initiated by self._init_analyse_timesteps
self.analyse_timesteps = None
+ # define array of timesteps at which to plot stuff
+ self.timesteps_to_plot = np.arange(0, self.number_of_timesteps, plot_timestep_every)
+ # we are now missing the last timestep.
+ self.timesteps_to_plot.append(self.number_of_timesteps)
# The number of subdomains are counted by self.init_meshes_and_markers()
self._number_of_subdomains = 0
# variable to check if self.set_parameters() has been called
@@ -338,6 +343,18 @@ class LDDsimulation(object):
df.info(colored("Start time stepping","green"))
# write initial values to file.
self.timestep_num = int(1)
+ # set up spacetime_errornorm dictionaries.
+ self.spacetime_errornorm = dict()
+ for subdom_ind, subdomain in self.subdomain.items():
+ self.spacetime_errornorm.update({subdom_ind: dict()})
+ for phase in subdomain.has_phases:
+ self.spacetime_errornorm[subdom_ind].update({phase: dict()})
+ self.spacetime_errornorm[subdom_ind][phase].update(
+ {# keys mean the norm in time. We use L2 in space
+ 'Linf': 0.0,
+ 'L1': 0.0,
+ 'L2': 0.0 }
+ )
# note that at this point of the code self.t = self.starttime as set in
# the beginning.
while self.timestep_num <= self.number_of_timesteps:
@@ -362,11 +379,12 @@ class LDDsimulation(object):
# for phase in subdomain.has_phases:
# print(f"calculated pressure of {phase}phase:")
# print(subdomain.pressure[phase].vector()[:])
- subdomain.write_solution_to_xdmf(#
- file = self.solution_file[subdom_ind], #
- time = self.t,#
- write_iter_for_fixed_time = False,
- )
+ if np.isin(self.timestep_num, self.timesteps_to_plot, assume_unique=True):
+ subdomain.write_solution_to_xdmf(#
+ file = self.solution_file[subdom_ind], #
+ time = self.t,#
+ write_iter_for_fixed_time = False,
+ )
# if we have an exact solution, write out the |u -uh|_L2 to see the
# absolute error.
if subdomain.pressure_exact is not None:
@@ -377,23 +395,11 @@ class LDDsimulation(object):
pa_exact.t = self.t
norm_pa_exact = df.norm(pa_exact, norm_type='L2', mesh=subdomain.mesh)
error_calculated = df.errornorm(pa_exact, subdomain.pressure[phase], 'L2', degree_rise=6)
- pressure_exact.update(
- {phase: df.interpolate(pa_exact, subdomain.function_space["pressure"][phase])}
- )
- absolute_difference = df.Function(subdomain.function_space["pressure"][phase])
- pressure_difference = pressure_exact[phase].vector()[:] - subdomain.pressure[phase].vector()[:]
- abs_diff_tmp = np.fabs(pressure_difference)
- absolute_difference.vector()[:] = abs_diff_tmp
- if phase == "wetting":
- data_set_name="abs_diff_w"
- else:
- data_set_name="abs_diff_nw"
- self.write_function_to_xdmf(#
- function=absolute_difference,
- file=self.solution_file[subdom_ind], #
- time=self.t,#
- data_set_name=data_set_name
- )
+ self.spacetime_errornorm[subdom_ind][phase]['L1'] += self.timestep_size*self.spacetime_errornorm[subdom_ind][phase]['L1']
+ self.spacetime_errornorm[subdom_ind][phase]['L2'] += self.timestep_size*self.spacetime_errornorm[subdom_ind][phase]['L2']**2
+ if error_calculated > self.spacetime_errornorm[subdom_ind][phase]['Linf']:
+ self.spacetime_errornorm[subdom_ind][phase]['Linf'] = error_calculated
+
if norm_pa_exact > self.tol:
relative_L2_errornorm.update(
{phase: error_calculated/norm_pa_exact}
@@ -402,15 +408,36 @@ class LDDsimulation(object):
relative_L2_errornorm.update(
{phase: error_calculated}
)
- errornorm_filename = self.output_dir+self.output_filename_parameter_part[subdom_ind]+\
- "_L2_errornorms_over_time" +".csv"
- self.write_errornorms_to_csv(
- filename = errornorm_filename, #
- subdomain_index = subdom_ind,
- errors = relative_L2_errornorm,
- )
+
+ if np.isin(self.timestep_num, self.timesteps_to_plot, assume_unique=True):
+ errornorm_filename = self.output_dir+self.output_filename_parameter_part[subdom_ind]+\
+ "_L2_errornorms_over_time" +".csv"
+ self.write_errornorms_to_csv(
+ filename = errornorm_filename, #
+ subdomain_index = subdom_ind,
+ errors = relative_L2_errornorm,
+ )
+ if np.isin(self.timestep_num, self.timesteps_to_plot, assume_unique=True):
+ pressure_exact.update(
+ {phase: df.interpolate(pa_exact, subdomain.function_space["pressure"][phase])}
+ )
+ absolute_difference = df.Function(subdomain.function_space["pressure"][phase])
+ pressure_difference = pressure_exact[phase].vector()[:] - subdomain.pressure[phase].vector()[:]
+ abs_diff_tmp = np.fabs(pressure_difference)
+ absolute_difference.vector()[:] = abs_diff_tmp
+ if phase == "wetting":
+ data_set_name="abs_diff_w"
+ else:
+ data_set_name="abs_diff_nw"
+ self.write_function_to_xdmf(#
+ function=absolute_difference,
+ file=self.solution_file[subdom_ind], #
+ time=self.t,#
+ data_set_name=data_set_name
+ )
self.timestep_num += 1
+ self.spacetime_errornorm[subdom_ind][phase]['L2'] = np.sqrt(self.spacetime_errornorm[subdom_ind][phase]['L2'])
print("LDD simulation use case: {}".format(self.use_case))
df.info("Finished calculation")
@@ -773,8 +800,7 @@ class LDDsimulation(object):
for subdom_ind, subdomain in self.subdomain.items():
t = copy.copy(self.starttime)
file = self.solution_file[subdom_ind]
-
- for timestep in range(self.number_of_timesteps+1):
+ for timestep in self.timesteps_to_plot:
exact_pressure = dict()
for phase in subdomain.has_phases:
pa_exact = subdomain.pressure_exact[phase]