diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index f9e7794580c0ab461dce1147d12b5119c56115cb..e9d60574be1fd1fb5d4e4996e8844cd4d4c85ab1 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -67,8 +67,8 @@ class LDDsimulation(object):
# df.parameters['std_out_all_processes'] = False
- # df.set_log_level(df.LogLevel.WARNING)
- df.set_log_level(df.LogLevel.INFO)
+ df.set_log_level(df.LogLevel.WARNING)
+ # df.set_log_level(df.LogLevel.INFO)
# df.set_log_level(df.LogLevel.DEBUG)
# CRITICAL = 50, // errors that may lead to data corruption and suchlike
@@ -155,19 +155,19 @@ class LDDsimulation(object):
# subdomain. This gets populated by self._init_solution_files()
self.solution_file = None
########## SOLVER DEFINITION AND PARAMETERS ########
- print("\nLinear Algebra Backends:")
- df.list_linear_algebra_backends()
- # # print("\nLinear Solver Methods:")
- df.list_linear_solver_methods()
- print("\n")
- df.list_krylov_solver_methods()
- print("\n")
- # # print("\nPeconditioners for Krylov Solvers:")
- df.list_krylov_solver_preconditioners()
+ # print("\nLinear Algebra Backends:")
+ # df.list_linear_algebra_backends()
+ # # # print("\nLinear Solver Methods:")
+ # df.list_linear_solver_methods()
+ # print("\n")
+ # df.list_krylov_solver_methods()
+ # print("\n")
+ # # # print("\nPeconditioners for Krylov Solvers:")
+ # df.list_krylov_solver_preconditioners()
- df.info(df.parameters, True)
+ # df.info(df.parameters, True)
self.solver_type_is_Krylov = True
### Define the linear solver to be used.
@@ -240,6 +240,8 @@ class LDDsimulation(object):
self.use_case = use_case
if include_gravity:
self.use_case = self.use_case + "-with-gravity"
+ if self.mesh_study:
+ self.use_case = "mesh study (meshres"+ "{}".format(mesh_resolution) + ") for case: " + self.use_case
self.output_dir = output_dir
self.subdomain_def_points = subdomain_def_points
self.isRichards = isRichards
@@ -259,10 +261,16 @@ class LDDsimulation(object):
self.number_of_timesteps = number_of_timesteps
self.number_of_timesteps_to_analyse = number_of_timesteps_to_analyse
# define array of timesteps at which to plot stuff
- self.timesteps_to_plot = np.arange(0, self.number_of_timesteps, plot_timestep_every)
+ self.timestep_numbers_to_plot = np.arange(0, self.number_of_timesteps, plot_timestep_every)
# we are now missing the last timestep.
- np.append(self.timesteps_to_plot, self.number_of_timesteps)
+ self.timestep_numbers_to_plot = np.append(self.timestep_numbers_to_plot, self.number_of_timesteps)
+ hlp.print_once("The following timesteps will be used for plotting:")
+ # print(self.timestep_numbers_to_plot)
self.timestep_size = timestep_size
+ # stores the timesteps tn that will be used to write out the solution.
+ # these are the timesteps corresponding to the timestep numbers in
+ # self.timestep_numbers_to_plot
+ self.timesteps_to_plot = self.starttime + self.timestep_size*self.timestep_numbers_to_plot
self.sources = sources
self.initial_conditions = initial_conditions
self.dirichletBC_expression_strings = dirichletBC_expression_strings
@@ -357,8 +365,9 @@ class LDDsimulation(object):
print(f"entering timestep ",
"{}".format(self.timestep_num),
"at time t = "+ "{number:.{digits}f}".format(number = self.t, digits = 4))
- # 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):
+
+ # check if the solver should be analised or not.
+ if np.isin(self.timestep_num, analyse_solver_at_timesteps, assume_unique=True):
print("analyising timestep ...")
self.LDDsolver(debug=self.debug_solver,
analyse_timestep=True,
@@ -369,71 +378,19 @@ class LDDsimulation(object):
analyse_condition=False)
self.t += self.timestep_size
- # write solutions to files.
- for subdom_ind in self.isRichards.keys():
- subdomain = self.subdomain[subdom_ind]
- # for phase in subdomain.has_phases:
- # print(f"calculated pressure of {phase}phase:")
- # print(subdomain.pressure[phase].vector()[:])
- if np.isin(self.timestep_num, self.timesteps_to_plot, assume_unique=True) and self.write2file['solutions']:
- 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:
- relative_L2_errornorm = dict()
- pressure_exact = dict()
- for phase in subdomain.has_phases:
- pa_exact = subdomain.pressure_exact[phase]
- 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)
- self.output[subdom_ind]['errornorm'][phase]['L1'] += self.timestep_size*error_calculated
- self.output[subdom_ind]['errornorm'][phase]['L2'] += self.timestep_size*error_calculated**2
- if error_calculated > self.output[subdom_ind]['errornorm'][phase]['Linf']:
- self.output[subdom_ind]['errornorm'][phase]['Linf'] = error_calculated
-
- if norm_pa_exact > self.tol:
- relative_L2_errornorm.update(
- {phase: error_calculated/norm_pa_exact}
- )
- else:
- relative_L2_errornorm.update(
- {phase: error_calculated}
- )
-
- 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) and self.write2file['absolute_differences']:
- 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
- )
+ # calculate spacetime errornorms for self.output.
+ space_errornorm = self._spacetime_errornorm_calc_step()
+ # write solutions to files. It is paramount, that self.timestep_num
+ # be only updated after calling self.write_data_to_disc()
+ self._write_data_to_disc(space_errornorm)
self.timestep_num += 1
- self.output[subdom_ind]['errornorm'][phase]['L2'] = np.sqrt(self.output[subdom_ind]['errornorm'][phase]['L2'])
+
+ # the L2 space_time errornorm needs to be square rooted.
+ for subdom_ind in self.isRichards.keys():
+ subdomain = self.subdomain[subdom_ind]
+ for phase in subdomain.has_phases:
+ self.output[subdom_ind]['errornorm'][phase]['L2'] = np.sqrt(self.output[subdom_ind]['errornorm'][phase]['L2'])
print("LDD simulation use case: {}".format(self.use_case))
df.info("Finished calculation")
@@ -453,6 +410,7 @@ class LDDsimulation(object):
# df.info(colored("finished post processing calculations \nAll right. I'm Done.", "green"))
return self.output
+
def LDDsolver(self, time: float = None, #
debug: bool = False, #
analyse_timestep: bool = False,
@@ -554,7 +512,7 @@ class LDDsimulation(object):
# of the subdomain for communication.
subdomain.write_pressure_to_interfaces()
- if analyse_timestep and self.write_to_file['L_iterations_per_timestep']:
+ if analyse_timestep and self.write2file['L_iterations_per_timestep']:
subdomain.write_solution_to_xdmf(#
file = solution_over_iteration_within_timestep[sd_index], #
time = time,#
@@ -567,17 +525,10 @@ class LDDsimulation(object):
)
# end loop over subdomains
##### stopping criterion for the solver.
- # only check if error criterion has been met after at least one iteration.
- #if iteration >= 1:
- 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}")
-
total_subsequent_error = np.amax(max_subsequent_error)
if debug:
- print(f"time = {time} (at timestep {self.timestep_num}):: max subsequent error of both subdomain in iteration {iteration} = ", total_subsequent_error)
+ print(f"the error criterion is tol = {error_criterion}")
+ print(f"time = {time} (at timestep {self.timestep_num}):: max subsequent error of both subdomains in iteration {iteration} is: ", total_subsequent_error)
if total_subsequent_error < error_criterion:
# if debug:
print(f"all phases on all subdomains reached a subsequent error = {total_subsequent_error} < tol = {error_criterion} after {iteration} iterations.")
@@ -744,6 +695,130 @@ class LDDsimulation(object):
return condition_number
+ def _spacetime_errornorm_calc_step(self)-> tp.Dict[int, tp.Dict[str, float]]:
+ """one step in integrating the spacetime errornorm
+
+ OUTPUT
+ error_out tp.Dict[int, tp.Dict[str, float]] dictionary containing
+ the L2-errornorms in space for the timestep.
+ this gets passed on to self.write_data_to_disc()
+ """
+
+ error_out = dict()
+ for subdom_ind in self.isRichards.keys():
+ error_out.update({subdom_ind: dict()})
+ subdomain = self.subdomain[subdom_ind]
+ # if we have an exact solution, write out several errors and
+ # errornorms to be used later in paraview or latex plots.
+ if subdomain.pressure_exact is not None:
+ for phase in subdomain.has_phases:
+ pa_exact = subdomain.pressure_exact[phase]
+ pa_exact.t = self.t
+ error_calculated = df.errornorm(pa_exact, subdomain.pressure[phase], 'L2', degree_rise=6)
+ self.output[subdom_ind]['errornorm'][phase]['L1'] += self.timestep_size*error_calculated
+ self.output[subdom_ind]['errornorm'][phase]['L2'] += self.timestep_size*error_calculated**2
+ if error_calculated > self.output[subdom_ind]['errornorm'][phase]['Linf']:
+ self.output[subdom_ind]['errornorm'][phase]['Linf'] = error_calculated
+
+ # if we are at a timestep at which to write shit out,
+ # calculate the relative errornorm
+ if np.isin(self.timestep_num, self.timestep_numbers_to_plot, assume_unique=True):
+ norm_pa_exact = df.norm(pa_exact, norm_type='L2', mesh=subdomain.mesh)
+ if norm_pa_exact > self.tol:
+ error_out[subdom_ind].update(
+ {phase: error_calculated/norm_pa_exact}
+ )
+ else:
+ error_out[subdom_ind].update(
+ {phase: error_calculated}
+ )
+ else:
+ error_out[subdom_ind].update({phase: None})
+ else:
+ # error_out.update({subdom_ind: None})
+ error_outupdate({subdom_ind: None})
+ return error_out
+
+
+ def _write_data_to_disc(self, space_errornorm: tp.Dict[int, tp.Dict[str, float]]):
+ """function that takes care of writing out different data to disc after
+ a timestep done by the self.run() function. What is written out at which
+ timesteps, is determined by by self.write2file aswell as
+ self.timestep_numbers_to_plot.
+
+ INPUT
+ space_errornorm tp.Dict[int, tp.Dict[str, float]]) dictionary of
+ L2 space errornorms calculated by
+ self._spacetime_errornorm_calc_step()
+ """
+ # we only write data of timesteps to disc, that are stored in
+ # self.timestep_numbers_to_plot.
+ if np.isin(self.timestep_num, self.timestep_numbers_to_plot, assume_unique=True):
+ # We need to determin the corresponding
+ # timestep value in self.timesteps_to_plot to write data out at the
+ # same values for timestep values tn. This is necessary, because we
+ # used self.timesteps_to_plot to plot the exact solutions.
+ # Counting up the time via self.t += self.timestep_size in the timesteping
+ # loop, yields to slightly other values for t due to rounding errors. sigh.
+ plot_index = np.argwhere(self.timestep_numbers_to_plot==self.timestep_num)[0][0]
+ plot_t = self.timesteps_to_plot[plot_index]
+
+ for subdom_ind in self.isRichards.keys():
+ subdomain = self.subdomain[subdom_ind]
+ # write out solution, if self.write2file['solutions'] flag says so.
+ if self.write2file['solutions']:
+ subdomain.write_solution_to_xdmf(#
+ file = self.solution_file[subdom_ind], #
+ time = plot_t,#
+ write_iter_for_fixed_time = False,
+ )
+ # if we have an exact solution, write out several errors and
+ # errornorms to be used later in paraview or latex plots.
+ if subdomain.pressure_exact is not None:
+ # write out the absolute differences for convenience in paraview
+ if self.write2file['absolute_differences']:
+ pressure_exact = dict()
+ for phase in subdomain.has_phases:
+ pa_exact = subdomain.pressure_exact[phase]
+ pa_exact.t = plot_t
+ 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=plot_t,#
+ data_set_name=data_set_name
+ )
+
+ if self.write2file['space_errornorms']:
+ relative_L2_errornorm = space_errornorm[subdom_ind]
+ 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,
+ time=plot_t,
+ timestep_num=self.timestep_num
+ )
+ else:
+ pass
+ pass
+ else:
+ # in this case we are at a timestep that is not bound to be written
+ # to disc.
+ pass
+
+
def _init_solution_files(self):
""" set up solution files for saving the solution of the LDD scheme for
each subdomain.
@@ -795,19 +870,19 @@ class LDDsimulation(object):
evaluate the exact solution of the simulation (in case there is one) at
all timesteps and write it to the solution file.
"""
+ # print(f" solution will be printed at times t = {self.timesteps_to_plot}\n")
for subdom_ind, subdomain in self.subdomain.items():
- t = copy.copy(self.starttime)
file = self.solution_file[subdom_ind]
for timestep in self.timesteps_to_plot:
exact_pressure = dict()
for phase in subdomain.has_phases:
pa_exact = subdomain.pressure_exact[phase]
- pa_exact.t = t
+ pa_exact.t = timestep
exact_pressure.update(
{phase: df.interpolate(pa_exact, subdomain.function_space["pressure"][phase])}
)
exact_pressure[phase].rename("exact_pressure_"+"{}".format(phase), "exact_pressure_"+"{}".format(phase))
- file.write(exact_pressure[phase], t)
+ file.write(exact_pressure[phase], timestep)
exact_capillary_pressure = df.Function(subdomain.function_space["pressure"]['wetting'])
if subdomain.isRichards:
@@ -816,8 +891,7 @@ class LDDsimulation(object):
pc_temp = exact_pressure["nonwetting"].vector().get_local() - exact_pressure["wetting"].vector().get_local()
exact_capillary_pressure.vector().set_local(pc_temp)
exact_capillary_pressure.rename("pc_exact", "pc_exact")
- file.write(exact_capillary_pressure, t)
- t += self.timestep_size
+ file.write(exact_capillary_pressure, timestep)
def write_subsequent_errors_to_csv(self,#
@@ -885,28 +959,24 @@ class LDDsimulation(object):
filename: str, #
subdomain_index: int,#
errors: tp.Dict[str, float],#
+ time: float,
+ timestep_num: int
) -> None:
""" write subsequent errors to csv file for plotting with pgfplots"""
- # in case no time is given, use the time of the class.
- subdomain = self.subdomain[subdomain_index]
- time = self.t
- timestep = self.timestep_num
-
subdomain = self.subdomain[subdomain_index]
- # iteration = subdomain.iteration_number
if subdomain.isRichards:
- self.errors = pd.DataFrame({ "time": self.t,
- "timestep": self.timestep_num,
+ self.errors = pd.DataFrame({ "time": time,
+ "timestep": timestep_num,
#"th. initial mass": df.assemble(rho*df.dx), \
#"energy": self.total_energy(), \
- "wetting": errors["wetting"]}, index=[timestep])
+ "wetting": errors["wetting"]}, index=[timestep_num])
else:
- self.errors = pd.DataFrame({ "time": self.t,
- "timestep": self.timestep_num,
+ self.errors = pd.DataFrame({ "time": time,
+ "timestep": timestep_num,
"wetting": errors["wetting"],
- "nonwetting": errors["nonwetting"]}, index=[timestep])
+ "nonwetting": errors["nonwetting"]}, index=[timestep_num])
# if self._rank == 0:
# check if file exists
if os.path.isfile(filename) == True:
@@ -1144,13 +1214,13 @@ class LDDsimulation(object):
self.subdomain[subdom_num].linear_solver = df.KrylovSolver(self.solver, self.preconditioner)
else:
self.subdomain[subdom_num].linear_solver = df.LUSolver()
- print("solver method:")
- for parameter, value in self.subdomain[subdom_num].linear_solver.parameters.items():
- print(f"\'{parameter}\': {value}")
+ # print("solver method:")
+ # for parameter, value in self.subdomain[subdom_num].linear_solver.parameters.items():
+ # print(f"\'{parameter}\': {value}")
# assign solver parameters set at the beginning of the LDDsimulation class.
solver_param = self.subdomain[subdom_num].linear_solver.parameters
- df.info(solver_param, True)
+ # df.info(solver_param, True)
for parameter_name in self.solver_parameters.keys():
solver_param[parameter_name] = self.solver_parameters[parameter_name]
@@ -1164,7 +1234,7 @@ class LDDsimulation(object):
filepath = self.output_dir+f"outer_boundary_marker_subdomain{subdom_num}.pvd"
df.File(filepath) << self.subdomain[subdom_num].outer_boundary_marker
- print(self.subdomain[subdom_num])
+ # print(self.subdomain[subdom_num])
def _get_interfaces(self, subdomain_index: int, #
@@ -1454,34 +1524,17 @@ class LDDsimulation(object):
def _init_analyse_timesteps(self):
- """ determin timesteps to anayise
-
- by setting
- """
- # We want to write out csv files containing subsequent L2 errors of iterations for
- # number_of_timesteps_to_analyse many timesteps. We want to do this with
- # analyse_timesteps = np.linspace(starttime, number_of_timesteps, number_of_timesteps_to_analyse)
- # see below. Therefor, we need to check if (number_of_timesteps_to_analyse-1) is
- # a divisor of number_of_timesteps,
+ """ determin timesteps to anayise """
if self.number_of_timesteps_to_analyse == 0:
self.analyse_timesteps = None
- elif self.number_of_timesteps_to_analyse == 1:
- self.analyse_timesteps = [1]
else:
- if not self.number_of_timesteps%(self.number_of_timesteps_to_analyse-1) == 0:
- # if number_of_timesteps%(number_of_timesteps_to_analyse-1) is not 0, we
- # round the division error to the nearest integer and redefine number_of_timesteps
- print(f"since number_of_timesteps/number_of_timesteps_to_analyse = {self.number_of_timesteps/self.number_of_timesteps_to_analyse}",
- "is not an integer,\n")
- new_analysing_interval = int(round(self.number_of_timesteps/(self.number_of_timesteps_to_analyse-1)))
- self.number_of_timesteps = new_analysing_interval*(self.number_of_timesteps_to_analyse)
- print(f"I'm resetting number_of_timesteps = {self.number_of_timesteps}\n")
self.analyse_timesteps = np.linspace(1, self.number_of_timesteps, self.number_of_timesteps_to_analyse, dtype=int)
- print(f"the following timesteps will be analysed: {self.analyse_timesteps}")
- self.endtime = self.number_of_timesteps*self.timestep_size
+
+ self.endtime = self.starttime + self.number_of_timesteps*self.timestep_size
# time_interval = [starttime, Tmax]
- print(f"The simulation interval is [{self.starttime}, {self.endtime}]")
+ print(f"The simulation interval is [{self.starttime}, {self.endtime}]\n")
+ print(f"the following timesteps will be analysed: {self.analyse_timesteps}\n")
def _init_simulation_output(self):