From f7fa571e3d967d7967ed21ab13141dadbfb6d98f Mon Sep 17 00:00:00 2001 From: David <forenkram@gmx.de> Date: Tue, 30 Jun 2020 15:16:09 +0200 Subject: [PATCH] make it possible to give initial iteration a timestep number --- LDDsimulation/LDDsimulation.py | 95 +++++++++++++++---- LDDsimulation/helpers.py | 1 + .../TP-R-layered_soil-all-params-one.py | 8 +- ...soil-g-but-same-perm-coarse-dt-longterm.py | 8 +- .../TP-R-layered_soil-g-but-same-perm.py | 24 +++-- .../layered_soil/TP-R-layered_soil.py | 6 +- .../mesh_studies/TP-R-2-patch-mesh-study.py | 9 +- 7 files changed, 111 insertions(+), 40 deletions(-) diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py index ef8e8ef..78921ff 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 eff9e51..7f3696b 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 2db3def..421a2af 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 a98285d..361e44c 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 013e15b..40c3700 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 babbf0e..880eb3a 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 ea48fc1..73de79c 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"]) -- GitLab