From 74472cc04eef918ef4e7ab0e3b0f16ff38827859 Mon Sep 17 00:00:00 2001
From: David Seus <david.seus@nians.uni-stuttgart.de>
Date: Sun, 8 Sep 2019 18:11:31 +0200
Subject: [PATCH] get meshstudy script to work
---
LDDsimulation/LDDsimulation.py | 35 ++++++---
...TP-TP-2-patch-pure-dd-convergence-study.py | 76 +++++++++++++------
2 files changed, 77 insertions(+), 34 deletions(-)
diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index 9260e68..d3c8c43 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -35,7 +35,7 @@ class LDDsimulation(object):
# maximal number of L-iterations that the LDD solver uses.
max_iter_num: int = 1000,
FEM_Lagrange_degree: int = 1,
- plot_timestep_every: int = 1):
+ mesh_study: bool = False):
hlp.print_once("\n ###############################################\n",
"############# Begin Simulation ################\n",
"###############################################\n")
@@ -97,6 +97,10 @@ class LDDsimulation(object):
self._max_iter_num = max_iter_num
# directory to write files output to.
self.output_dir = None
+ # toggle wether or not we are doing a mesh study. influences the naming
+ # of the output directory
+ self.mesh_study = mesh_study
+
# list of subdomain meshes. Is populated by self.init_meshes_and_markers()
self.mesh_subdomain = []
# global marker for subdomains and interfaces. Is set by self.init_meshes_and_markers()
@@ -116,10 +120,6 @@ 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
@@ -234,6 +234,7 @@ class LDDsimulation(object):
densities: tp.Dict[int, tp.Dict[str, float]] = None,#
include_gravity: bool = False,#
gravity_acceleration: float = 9.81,
+ plot_timestep_every: int = 1,
)-> None:
@@ -260,6 +261,10 @@ class LDDsimulation(object):
# 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.timesteps_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_size = timestep_size
self.sources = sources
self.initial_conditions = initial_conditions
@@ -344,7 +349,7 @@ class LDDsimulation(object):
df.info(colored("Start time stepping","green"))
# write initial values to file.
self.timestep_num = int(1)
-
+
# 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:
@@ -365,7 +370,8 @@ class LDDsimulation(object):
self.t += self.timestep_size
# write solutions to files.
- for subdom_ind, subdomain in self.subdomain.items():
+ 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()[:])
@@ -385,8 +391,8 @@ 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)
- self.output[subdom_ind]['errornorm'][phase]['L1'] += self.timestep_size*self.output[subdom_ind]['errornorm'][phase]['L1']
- self.output[subdom_ind]['errornorm'][phase]['L2'] += self.timestep_size*self.output[subdom_ind]['errornorm'][phase]['L2']**2
+ 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
@@ -1440,7 +1446,12 @@ class LDDsimulation(object):
gravity_string = ""
if self.include_gravity:
gravity_string = "_with_gravity"
- dirname = "mesh_res" + "{}_".format(self.mesh_resolution)\
+ # if we want to investigate the EOC and mesh dependence, we do not want
+ # to put all data files with different grids in different directories.
+ if self.mesh_study:
+ dirname = gravity_string + "/"
+ else:
+ dirname = "mesh_res" + "{}_".format(self.mesh_resolution)\
+ "dt" + "{}".format(self.timestep_size)+ gravity_string + "/"
self.output_dir += dirname
@@ -1478,8 +1489,8 @@ class LDDsimulation(object):
self.endtime = self.number_of_timesteps*self.timestep_size
# time_interval = [starttime, Tmax]
print(f"The simulation interval is [{self.starttime}, {self.endtime}]")
-
-
+
+
def _init_simulation_output(self):
"""set up dictionary for simulation output"""
self.output = dict()
diff --git a/TP-TP-2-patch-pure-dd-avoid-interface-at-origin/mesh_study_convergence/TP-TP-2-patch-pure-dd-convergence-study.py b/TP-TP-2-patch-pure-dd-avoid-interface-at-origin/mesh_study_convergence/TP-TP-2-patch-pure-dd-convergence-study.py
index a4d2130..01fc996 100755
--- a/TP-TP-2-patch-pure-dd-avoid-interface-at-origin/mesh_study_convergence/TP-TP-2-patch-pure-dd-convergence-study.py
+++ b/TP-TP-2-patch-pure-dd-avoid-interface-at-origin/mesh_study_convergence/TP-TP-2-patch-pure-dd-convergence-study.py
@@ -9,6 +9,8 @@ import LDDsimulation as ldd
import functools as ft
import helpers as hlp
import datetime
+import os
+import pandas as pd
date = datetime.datetime.now()
datestr = date.strftime("%Y-%m-%d")
@@ -21,12 +23,13 @@ use_case = "TP-TP-2-patch-pure-dd"
solver_tol = 6E-7
max_iter_num = 1000
FEM_Lagrange_degree = 1
+mesh_study = True
############ GRID #######################
# mesh_resolution = 20
timestep_size = 0.0001
-number_of_timesteps = 1000
-plot_timestep_every = 10
+number_of_timesteps = 100
+plot_timestep_every = 1
# decide how many timesteps you want analysed. Analysed means, that we write out
# subsequent errors of the L-iteration within the timestep.
number_of_timesteps_to_analyse = 0
@@ -39,8 +42,8 @@ lambda_w = 40
lambda_nw = 40
include_gravity = False
-debugflag = False
-analyse_condition = True
+debugflag = True
+analyse_condition = False
output_string = "./output/{}-{}_timesteps{}_P{}-solver_tol{}".format(datestr, use_case, number_of_timesteps, FEM_Lagrange_degree, solver_tol)
@@ -510,14 +513,24 @@ for subdomain in isRichards.keys():
#
# sa
-write_to_file = {
- 'meshes_and_markers': True,
- 'L_iterations': True,
- 'solutions': False,
- 'absolute_differences': False,
- 'condition_numbers': False,
- 'subsequent_errors': False
-}
+if mesh_study:
+ write_to_file = {
+ 'meshes_and_markers': True,
+ 'L_iterations': True,
+ 'solutions': False,
+ 'absolute_differences': False,
+ 'condition_numbers': False,
+ 'subsequent_errors': False
+ }
+else:
+ write_to_file = {
+ 'meshes_and_markers': True,
+ 'L_iterations': True,
+ 'solutions': True,
+ 'absolute_differences': True,
+ 'condition_numbers': True,
+ 'subsequent_errors': True
+ }
# initialise LDD simulation class
@@ -527,10 +540,10 @@ simulation = ldd.LDDsimulation(
debug=debugflag,
max_iter_num=max_iter_num,
FEM_Lagrange_degree=FEM_Lagrange_degree,
- plot_timestep_every=plot_timestep_every
+ mesh_study=mesh_study
)
-resolutions = [5, 10, 15] #, 20, 25, 30, 35, 40, 45, 50]
+resolutions = [5, 10, 15, 20] #, 20, 25, 30, 35, 40, 45, 50]
for mesh_resolution in resolutions:
simulation.set_parameters(use_case=use_case,
@@ -550,6 +563,7 @@ for mesh_resolution in resolutions:
starttime=starttime,
number_of_timesteps=number_of_timesteps,
number_of_timesteps_to_analyse=number_of_timesteps_to_analyse,
+ plot_timestep_every=plot_timestep_every,
timestep_size=timestep_size,
sources=source_expression,
initial_conditions=initial_condition,
@@ -561,13 +575,31 @@ for mesh_resolution in resolutions:
)
simulation.initialise()
+ output_dir = simulation.output_dir
# simulation.write_exact_solution_to_xdmf()
output = simulation.run(analyse_condition=analyse_condition)
- for subdomain_index in output.keys():
- mesh_h = output['mesh_size']
- for phase, different_errornorms in output[subdomain_index]['errornorm'].items():
- for errortype, errornorm in different_errornorms.items():
- eoc_filename = "{}_{}-errornorm-for-{}".format(output_string, errortype, phase)
- eocfile = open("eoc_filename", "a")
- eocfile.write( str(mesh_h) + " " + str(errornorm) + "\n" )
- eocfile.close()
+ for subdomain_index, subdomain_output in output.items():
+ mesh_h = subdomain_output['mesh_size']
+ for phase, different_errornorms in subdomain_output['errornorm'].items():
+ filename = output_dir + "subdomain{}-space-time-errornorm-{}-phase.csv".format(subdomain_index, phase)
+ # for errortype, errornorm in different_errornorms.items():
+
+ # eocfile = open("eoc_filename", "a")
+ # eocfile.write( str(mesh_h) + " " + str(errornorm) + "\n" )
+ # eocfile.close()
+ # if subdomain.isRichards:mesh_h
+ data_dict = {
+ 'mesh_parameter': mesh_resolution,
+ 'mesh_h': mesh_h,
+ }
+ for error_type, errornorms in different_errornorms.items():
+ data_dict.update(
+ {error_type: errornorms}
+ )
+ errors = pd.DataFrame(data_dict, index=[mesh_resolution])
+ # check if file exists
+ if os.path.isfile(filename) == True:
+ with open(filename, 'a') as f:
+ errors.to_csv(f, header=False, sep='\t', encoding='utf-8', index=False)
+ else:
+ errors.to_csv(filename, sep='\t', encoding='utf-8', index=False)
--
GitLab