diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index ef8e8efd749fed7fc55d7415b0bc435fec1e2143..78921ffd83504dd8dc5d4745ca33edcb53d5a0b6 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -235,6 +235,7 @@ class LDDsimulation(object):
densities: tp.Dict[int, tp.Dict[str, float]] = None,#
include_gravity: bool = False,#
gravity_acceleration: float = 9.81,
+ starttime_timestep_number_shift: int = 0,
plot_timestep_every: int = 1,
)-> None:
@@ -260,20 +261,40 @@ class LDDsimulation(object):
self.saturation = saturation
# simulation start time and time variable
self.starttime = starttime
+ # optional argument specifying what number the initial timestep should
+ # have. Default is 0, but if a simulation needs to be restarted from
+ # a later timestep, it is convenient for merging data to shift the
+ # the numbering of the timestep by the sepcified amount.
+ self.starttime_timestep_number_shift = starttime_timestep_number_shift
# total number of timesteps.
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.timestep_numbers_to_plot = np.arange(0, self.number_of_timesteps, plot_timestep_every)
+ # define array of timesteps numbers for which to plot stuff. This
+ # does not yet contain the shift self.starttime_timestep_number_shift,
+ # because the nonshifted array is needed to calculate
+ # self.timesteps_to_plot below.
+ self.timestep_numbers_to_plot = np.arange(
+ 0,
+ self.number_of_timesteps,
+ plot_timestep_every
+ )
# we are now missing the last timestep.
- self.timestep_numbers_to_plot = np.append(self.timestep_numbers_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.timesteps_to_plot = self.starttime \
+ + self.timestep_size*self.timestep_numbers_to_plot
+ # now the numbering of the timesteps to plot can safely be shifted by
+ # the optional argument self.starttime_timestep_number_shift
+ self.timestep_numbers_to_plot = self.timestep_numbers_to_plot \
+ + self.starttime_timestep_number_shift
self.sources = sources
self.initial_conditions = initial_conditions
self.dirichletBC_expression_strings = dirichletBC_expression_strings
@@ -358,12 +379,14 @@ class LDDsimulation(object):
df.info(colored("Start time stepping","green"))
# write initial values to file.
- self.timestep_num = int(1)
+ self.timestep_num = int(1) + self.starttime_timestep_number_shift
self.t = self.starttime
# 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:
+ self.max_timestep_num = self.number_of_timesteps\
+ + self.starttime_timestep_number_shift
+ while self.timestep_num <= self.max_timestep_num:
print("LDD simulation use case: {}".format(self.use_case))
print(f"entering timestep ",
"{}".format(self.timestep_num),
@@ -1056,41 +1079,69 @@ class LDDsimulation(object):
if self.exact_solution:
exact_pressure = dict()
S = self.saturation[subdom_ind]
- saturation_w = df.Function(subdomain.function_space["pressure"]['wetting'])
- saturation_nw = df.Function(subdomain.function_space["pressure"]['wetting'])
+ saturation_w = df.Function(
+ subdomain.function_space["pressure"]['wetting']
+ )
+ saturation_nw = df.Function(
+ subdomain.function_space["pressure"]['wetting']
+ )
for phase in subdomain.has_phases:
f_expr = subdomain.source[phase]
f_expr.t = timestep
source.update(
- {phase: df.project(f_expr, subdomain.function_space["pressure"][phase])}
+ {phase: df.project(
+ f_expr, subdomain.function_space["pressure"][phase]
+ )
+ }
)
- source[phase].rename("source_"+"{}".format(phase), "source_"+"{}".format(phase))
+ source[phase].rename(
+ "source_"+"{}".format(phase),
+ "source_"+"{}".format(phase)
+ )
file.write(source[phase], timestep)
if self.exact_solution:
pa_exact = subdomain.pressure_exact[phase]
pa_exact.t = timestep
exact_pressure.update(
- {phase: df.project(pa_exact, subdomain.function_space["pressure"][phase])}
+ {phase: df.project(
+ pa_exact,
+ subdomain.function_space["pressure"][phase]
+ )}
+ )
+ exact_pressure[phase].rename(
+ "exact_pressure_"+"{}".format(phase),
+ "exact_pressure_"+"{}".format(phase)
)
- exact_pressure[phase].rename("exact_pressure_"+"{}".format(phase), "exact_pressure_"+"{}".format(phase))
file.write(exact_pressure[phase], timestep)
if self.exact_solution:
- exact_capillary_pressure = df.Function(subdomain.function_space["pressure"]['wetting'])
+ exact_capillary_pressure = df.Function(
+ subdomain.function_space["pressure"]['wetting']
+ )
if subdomain.isRichards:
- exact_capillary_pressure.assign(-exact_pressure['wetting'])
+ exact_capillary_pressure.assign(
+ -exact_pressure['wetting']
+ )
else:
- pc_temp = exact_pressure["nonwetting"].vector().get_local() - exact_pressure["wetting"].vector().get_local()
+ 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, timestep)
- saturation_w = df.project(S(exact_capillary_pressure), subdomain.function_space["pressure"]["wetting"])
+ saturation_w = df.project(
+ S(exact_capillary_pressure),
+ subdomain.function_space["pressure"]["wetting"]
+ )
# saturation_w.assign(Sat_w)
saturation_w.rename("Sw_exact", "Sw_exact")
file.write(saturation_w, timestep)
# S_nw = 1-S(exact_capillary_pressure).vector().get_local()
- saturation_nw = df.project(1-S(exact_capillary_pressure), subdomain.function_space["pressure"]["wetting"])
+ saturation_nw = df.project(
+ 1-S(exact_capillary_pressure),
+ subdomain.function_space["pressure"]["wetting"]
+ )
# saturation_nw.assign(S_nw)
saturation_nw.rename("Snw_exact", "Snw_exact")
file.write(saturation_nw, timestep)
@@ -1739,10 +1790,16 @@ class LDDsimulation(object):
if self.number_of_timesteps_to_analyse == 0:
self.analyse_timesteps = None
else:
- self.analyse_timesteps = np.linspace(1, self.number_of_timesteps, self.number_of_timesteps_to_analyse, dtype=int)
+ self.analyse_timesteps = np.linspace(
+ 1 + self.starttime_timestep_number_shift,
+ self.number_of_timesteps + self.starttime_timestep_number_shift,
+ self.number_of_timesteps_to_analyse,
+ dtype=int
+ )
- self.endtime = self.starttime + 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}]\n")
print(f"the following timesteps will be analysed: {self.analyse_timesteps}\n")
diff --git a/LDDsimulation/helpers.py b/LDDsimulation/helpers.py
index eff9e51347e84bde83d9b4b5092c83d976048751..7f3696b47096540efdb966293c79de9da8d3a193 100644
--- a/LDDsimulation/helpers.py
+++ b/LDDsimulation/helpers.py
@@ -61,6 +61,7 @@ def run_simulation(
intrinsic_permeability=parameter["intrinsic_permeability"],
saturation=parameter["sat_pressure_relationship"],
starttime=starttime, # parameter["starttime"],
+ starttime_timestep_number_shift=parameter["starttime_timestep_number_shift"],
number_of_timesteps=parameter["number_of_timesteps"],
number_of_timesteps_to_analyse=parameter["number_of_timesteps_to_analyse"],
plot_timestep_every=parameter["plot_timestep_every"],
diff --git a/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-all-params-one.py b/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-all-params-one.py
index 2db3def978f40226b2e3fa1c08f47d6681585688..421a2af49040c4ca411f937850c61ab5b75f354e 100755
--- a/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-all-params-one.py
+++ b/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-all-params-one.py
@@ -57,7 +57,7 @@ resolutions = {
# starttimes gives a list of starttimes to run the simulation from.
# The list is looped over and a simulation is run with t_0 as initial time
# for each element t_0 in starttimes.
-starttimes = [0.0]
+starttimes = {0: 0.0}
timestep_size = 0.001
number_of_timesteps = 5
@@ -325,7 +325,6 @@ f = open(thisfile, 'r')
print(f.read())
f.close()
-
# MAIN ########################################################################
if __name__ == '__main__':
# dictionary of simualation parameters to pass to the run function.
@@ -368,7 +367,10 @@ if __name__ == '__main__':
"write_to_file": write_to_file,
"analyse_condition": analyse_condition
}
- for starttime in starttimes:
+ for number_shift, starttime in starttimes.items():
+ simulation_parameter.update(
+ {"starttime_timestep_number_shift": number_shift}
+ )
for mesh_resolution, solver_tol in resolutions.items():
simulation_parameter.update({"solver_tol": solver_tol})
hlp.info(simulation_parameter["use_case"])
diff --git a/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-g-but-same-perm-coarse-dt-longterm.py b/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-g-but-same-perm-coarse-dt-longterm.py
index a98285d55cc2cdb43245eae7ea37116ffeb989fe..361e44c05c6243c5e1f44bac14f5cbddfc2d62e1 100755
--- a/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-g-but-same-perm-coarse-dt-longterm.py
+++ b/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-g-but-same-perm-coarse-dt-longterm.py
@@ -57,7 +57,7 @@ resolutions = {
# starttimes gives a list of starttimes to run the simulation from.
# The list is looped over and a simulation is run with t_0 as initial time
# for each element t_0 in starttimes.
-starttimes = [0.0]
+starttimes = {0: 0.0}
timestep_size = 0.01
number_of_timesteps = 400
@@ -340,7 +340,6 @@ f = open(thisfile, 'r')
print(f.read())
f.close()
-
# MAIN ########################################################################
if __name__ == '__main__':
# dictionary of simualation parameters to pass to the run function.
@@ -383,7 +382,10 @@ if __name__ == '__main__':
"write_to_file": write_to_file,
"analyse_condition": analyse_condition
}
- for starttime in starttimes:
+ for number_shift, starttime in starttimes.items():
+ simulation_parameter.update(
+ {"starttime_timestep_number_shift": number_shift}
+ )
for mesh_resolution, solver_tol in resolutions.items():
simulation_parameter.update({"solver_tol": solver_tol})
hlp.info(simulation_parameter["use_case"])
diff --git a/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-g-but-same-perm.py b/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-g-but-same-perm.py
index 013e15beca6f60b1f7c5713e1a6e5e81d4940dee..40c37006fb396446f3fbf081c9a582d60a3c532a 100755
--- a/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-g-but-same-perm.py
+++ b/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-g-but-same-perm.py
@@ -31,7 +31,7 @@ datestr = date.strftime("%Y-%m-%d")
# Name of the usecase that will be printed during simulation.
-use_case = "TP-R-layered-soil-realistic-same-perm"
+use_case = "TP-R-layered-soil-realistic-g-same-perm"
# The name of this very file. Needed for creating log output.
thisfile = "TP-R-layered_soil-g-but-same-perm.py"
@@ -54,12 +54,15 @@ resolutions = {
# 256: 1e-6,
}
-# starttimes gives a list of starttimes to run the simulation from.
-# The list is looped over and a simulation is run with t_0 as initial time
-# for each element t_0 in starttimes.
-starttimes = [0.0]
+# starttimes gives a dictionary of starttimes to run the simulation from.
+# The dictionary is looped over and a simulation is run for each item
+# timestep_num: t_0, t_0 being the initial time and timestep_num the number
+# that this initial timestep should be counted as.
+# Format:
+# timestep_num: initial_time
+starttimes = {650: 0.650}
timestep_size = 0.001
-number_of_timesteps = 800
+number_of_timesteps = 550
# LDD scheme parameters ######################################################
@@ -87,10 +90,10 @@ analyse_condition = False
# when number_of_timesteps is high, it might take a long time to write all
# timesteps to disk. Therefore, you can choose to only write data of every
# plot_timestep_every timestep to disk.
-plot_timestep_every = 4
+plot_timestep_every = 1
# Decide how many timesteps you want analysed. Analysed means, that
# subsequent errors of the L-iteration within the timestep are written out.
-number_of_timesteps_to_analyse = 5
+number_of_timesteps_to_analyse = 1
# fine grained control over data to be written to disk in the mesh study case
# as well as for a regular simuation for a fixed grid.
@@ -373,7 +376,10 @@ if __name__ == '__main__':
"write_to_file": write_to_file,
"analyse_condition": analyse_condition
}
- for starttime in starttimes:
+ for number_shift, starttime in starttimes.items():
+ simulation_parameter.update(
+ {"starttime_timestep_number_shift": number_shift}
+ )
for mesh_resolution, solver_tol in resolutions.items():
simulation_parameter.update({"solver_tol": solver_tol})
hlp.info(simulation_parameter["use_case"])
diff --git a/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil.py b/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil.py
index babbf0e2f639be1a2a1131943b33dcbcb59c8681..880eb3a6ef56a8c4e7c172c551468d685926f59f 100755
--- a/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil.py
+++ b/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil.py
@@ -337,7 +337,6 @@ f = open(thisfile, 'r')
print(f.read())
f.close()
-
# MAIN ########################################################################
if __name__ == '__main__':
# dictionary of simualation parameters to pass to the run function.
@@ -380,7 +379,10 @@ if __name__ == '__main__':
"write_to_file": write_to_file,
"analyse_condition": analyse_condition
}
- for starttime in starttimes:
+ for number_shift, starttime in starttimes.items():
+ simulation_parameter.update(
+ {"starttime_timestep_number_shift": number_shift}
+ )
for mesh_resolution, solver_tol in resolutions.items():
simulation_parameter.update({"solver_tol": solver_tol})
hlp.info(simulation_parameter["use_case"])
diff --git a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study.py
index ea48fc1c03542a282e2059bdd7b96aafbdaca244..73de79ca9cabfbfa75f9e54a65a0db7d400e179c 100755
--- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study.py
+++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study.py
@@ -56,7 +56,7 @@ resolutions = {
# starttimes gives a list of starttimes to run the simulation from.
# The list is looped over and a simulation is run with t_0 as initial time
# for each element t_0 in starttimes.
-starttimes = [0.0]
+starttimes = {0: 0.0}
timestep_size = 0.001
number_of_timesteps = 1000
@@ -170,7 +170,6 @@ L = {#
'nonwetting': Lnw2}
}
-
lambda_param = {#
# subdom_num : lambda parameter for the L-scheme
0 : {'wetting' :lambda_w,
@@ -267,7 +266,6 @@ f = open(thisfile, 'r')
print(f.read())
f.close()
-
# MAIN ########################################################################
if __name__ == '__main__':
# dictionary of simualation parameters to pass to the run function.
@@ -310,7 +308,10 @@ if __name__ == '__main__':
"write_to_file": write_to_file,
"analyse_condition": analyse_condition
}
- for starttime in starttimes:
+ for number_shift, starttime in starttimes.items():
+ simulation_parameter.update(
+ {"starttime_timestep_number_shift": number_shift}
+ )
for mesh_resolution, solver_tol in resolutions.items():
simulation_parameter.update({"solver_tol": solver_tol})
hlp.info(simulation_parameter["use_case"])