From d7b70a4ed06b359cc1bae47571a731c582a42492 Mon Sep 17 00:00:00 2001
From: David Seus <david.seus@ians.uni-stuttgart.de>
Date: Fri, 19 Apr 2019 16:50:17 +0200
Subject: [PATCH] write output for L2 errornorm with exact solution
---
LDDsimulation/LDDsimulation.py | 185 +++++++--
LDDsimulation/domainPatch.py | 16 +-
RR-2-patch-test-case/RR-2-patch-test.py | 381 ++----------------
.../plots/exact_solution_error_norms.tex | 47 +++
.../plots/subsequent_errors.tex | 114 +++++-
TP-TP-patch-test-case/TP-TP-2-patch-test.py | 45 ++-
6 files changed, 386 insertions(+), 402 deletions(-)
create mode 100644 RR-2-patch-test-case/plots/exact_solution_error_norms.tex
diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index 7e38bcd..74d5827 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -121,9 +121,12 @@ class LDDsimulation(object):
# dictionary to store the df.dirichletBC objects in. These are used by the
# solver to set the boundary conditions
self.outerBC = dict()
- # Flag to toggle between a case with exact solution and a purely numerical
- # simulation. This changes the behaviour of self.update_DirichletBC_dictionary()
- self.exact_solution_case = False
+ # Dictionary that holdes the exact solution. Gets set by self.set_parameters()
+ # but can remain None if we don't consider a case with an exact solution.
+ self.exact_solution = None
+ # Dictionary of df. Expressions holding for each subdomain and phase the
+ # dolfin expressions generated of self.exact_solution by self._init_initial_values()
+ self.exact_solution_expr = dict()
# dictionary holding the source expressions for each subdomain. These are
# saved to the respective subdomains by the self._init_subdomains() method.
self.sources = None#
@@ -178,7 +181,7 @@ class LDDsimulation(object):
initial_conditions: tp.Dict[int, tp.Dict[str, str]],#
dirichletBC_Expressions: tp.Dict[int, tp.Dict[str, str]],#
write2file: tp.Dict[str, bool],#
- exact_solution_case: bool = False, #
+ exact_solution: tp.Dict[int, tp.Dict[str, str]] = None, #
)-> None:
""" set parameters of an instance of class LDDsimulation"""
@@ -196,6 +199,7 @@ class LDDsimulation(object):
self.relative_permeability = relative_permeability
self.saturation = saturation
# simulation start time and time variable
+ self.t_initial = time_interval[0]
self.t = time_interval[0]
# end_time
self.T = time_interval[1]
@@ -203,7 +207,7 @@ class LDDsimulation(object):
self.sources = sources
self.initial_conditions = initial_conditions
self.dirichletBC_Expressions = dirichletBC_Expressions
- self.exact_solution_case = exact_solution_case
+ self.exact_solution = exact_solution
self.write2file = write2file
self._parameters_set = True
@@ -214,6 +218,7 @@ class LDDsimulation(object):
"""
if not self._parameters_set:
raise RuntimeError('Parameters note set!\nYou forgott to run self.set_parameters(**kwds)')
+ self._add_parameters_to_output_dirname()
print("initialise meshes and marker functions ...\n")
self._init_meshes_and_markers()
print("initialise interfaces ...\n")
@@ -225,9 +230,12 @@ class LDDsimulation(object):
# initial values get initialised here to be able to call the solver at different
# timesteps without having to call self.run().
self._init_initial_values()
+ # initialise exact solution expression dictionary if we have an exact
+ # solution case.
+ self._init_exact_solution_expression()
self._initialised = True
- def run(self, analyse_solver_at_times: tp.List[float] = None):
+ def run(self, analyse_solver_at_timesteps: tp.List[float] = None):
""" run the simulation
This method implements the time stepping, calls the acutal solvers and writes
@@ -235,7 +243,7 @@ class LDDsimulation(object):
want to understand this code, start reading here.
PARAMETERS
- analyse_solver_at_times #type tp.List[float] List of timesteps at which
+ analyse_solver_at_timesteps #type tp.List[float] List of timesteps at which
to produce output to
analyse the solver at
given timestep. This
@@ -256,15 +264,19 @@ class LDDsimulation(object):
df.info(colored("Start time stepping","green"))
# write initial values to file.
+ self.timestep_num = 0
while self.t < self.T:
- print(f"entering timestep t = self.t")
+ print(f"entering timestep ",
+ "{number:.{digits}f}".format(number = self.timestep_num, digits = 4),
+ f"t = {self.t}")
# check if the solver should beanalised or not.
- if analyse_solver_at_times is not None and np.isin(self.t, analyse_solver_at_times):
+ if analyse_solver_at_timesteps is not None and np.isin(self.timestep_num, analyse_solver_at_timesteps):
self.LDDsolver(debug = False, analyse_timestep = True)
else:
self.LDDsolver(debug = False, analyse_timestep = False)
self.t += self.timestep_size
+ self.timestep_num += 1
# write solutions to files.
for subdom_ind, subdomain in self.subdomain.items():
subdomain.write_solution_to_xdmf(#
@@ -272,6 +284,28 @@ class LDDsimulation(object):
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 self.exact_solution is not None:
+ L2errornorm = dict()
+ for phase in subdomain.has_phases:
+ #error = df.Function(subdomain.function_space[phase])
+ # error.assign(subdomain.pressure[phase] - subdomain.pressure_prev_iter[phase])
+ pa_exact = self.exact_solution_expr[subdom_ind][phase]
+ pa_exact.t = self.t
+ error_calculated = df.errornorm(pa_exact, subdomain.pressure[phase], 'L2', degree_rise=6)
+ L2errornorm.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 = L2errornorm,
+ )
+
+
df.info("Finished")
df.info("Elapsed time:"+str(timer.elapsed()[0]))
@@ -311,7 +345,9 @@ class LDDsimulation(object):
# create xdmf files for writing out iterations if analyse_timestep is
# given.
if analyse_timestep:
- filename = self.output_dir+"subdomain"+"{}_".format(ind)+ "iterations_at_time"+ "{}".format(time) +".xdmf"
+ filename = self.output_dir+self.output_filename_parameter_part[ind]\
+ + "solution_iterations_at_t"\
+ + "{number:.{digits}f}".format(number=time, digits=2) +".xdmf"
solution_over_iteration_within_timestep.update( #
{ind: SolutionFile(self._mpi_communicator, filename)}
)
@@ -376,26 +412,34 @@ class LDDsimulation(object):
subsequent_iter_error.update(
{phase: error_calculated}
)
- print(f"the subsequent error in iteration {iteration} is {subsequent_iter_error[phase]}")
+ if debug:
+ print(f"the subsequent error in iteration {iteration} is {subsequent_iter_error[phase]}")
if analyse_timestep:
- #
- subsequent_errors_over_iteration[sd_index].update(
- { phase: error })
- subsequent_error_filename = self.output_dir+"subdomain"+"{}_".format(sd_index)+\
- "subsequent_iteration_error_for_phase_"+"{}".format(phase)+"_at_time"+ "{}".format(time) +".csv"
- self.write_subsequent_errors(
+ # save the calculated L2 errors to csv file for plotting
+ # with pgfplots.
+ print(f"t = {self.t} the subsequent error in iteration {iteration} is {subsequent_iter_error[phase]}")
+ subsequent_error_filename = self.output_dir\
+ +self.output_filename_parameter_part[sd_index]\
+ +"subsequent_iteration_errors" +"_at_time"+\
+ "{number:.{digits}f}".format(number=time, digits=2) +".csv"
+ self.write_subsequent_errors_to_csv(
filename = subsequent_error_filename, #
subdomain_index = sd_index,
- subsequent_errors = subsequent_iter_error,
+ errors = subsequent_iter_error,
time = time)
+ # this saves the actual difference function between subsequent
+ # iterations to analyse then in paraview later.
+ subsequent_errors_over_iteration[sd_index].update(
+ { phase: error })
+
# both phases should be done before we mark the subdomain as
# finished.
total_subsequent_iter_error = max(subsequent_iter_error.values())
# check if error criterion has been met
- error_criterion = subdomain.mesh.hmin()**2*subdomain.timestep_size
+ error_criterion = min(1E-9, subdomain.mesh.hmin()**2*subdomain.timestep_size)
if debug:
print(f" total_subsequent_error in iter {iteration} = {total_subsequent_iter_error }\n")
print(f"the grid size on subdomain{sd_index} is h={subdomain.mesh.hmin()}")
@@ -440,6 +484,7 @@ class LDDsimulation(object):
iteration += 1
# end iteration while loop.
+
def Lsolver_step(self, subdomain_index: int,#
debug: bool = False) -> None:
""" L-scheme solver iteration step for an object of class subdomain
@@ -518,6 +563,7 @@ class LDDsimulation(object):
each subdomain.
"""
# set up solution file names for each subdomain.
+ self.output_filename_parameter_part = dict()
self.solution_file = dict()
for subdom_ind, subdomain in self.subdomain.items():
tau = subdomain.timestep_size
@@ -525,32 +571,36 @@ class LDDsimulation(object):
h = subdomain.mesh.hmin()
if subdomain.isRichards:
Lam = subdomain.lambda_param['wetting']
- name1 ="solution_subdomain" + "{}_".format(subdom_ind)
+ name1 ="subdomain" + "{}".format(subdom_ind)
# only show three digits of the mesh size.
- name2 = "dt" + "{}".format(tau)+"_hmin" + "{number:.{digits}f}".format(number=h, digits=3)
+ name2 = "_dt" + "{}".format(tau)+"_hmin" + "{number:.{digits}f}".format(number=h, digits=3)
name3 = "_lambda" + "{}".format(Lam) + "_Lw" + "{}".format(L["wetting"])
- filename = self.output_dir + name1 + name2 + name3 +".xdmf"
+ filename_param_part = name1 + name2 + name3
+ filename = self.output_dir + filename_param_part + "_solution" +".xdmf"
else:
Lam_w = subdomain.lambda_param['wetting']
Lam_nw = subdomain.lambda_param['nonwetting']
- filename = self.output_dir+"solution_subdomain" + "{}_".format(subdom_ind)\
- +"dt" + "{}".format(tau)\
+ filename_param_part = "subdomain" + "{}".format(subdom_ind)\
+ +"_dt" + "{}".format(tau)\
+"_hmin" + "{number:.{digits}f}".format(number=h, digits=3)\
+"_lambda_w" + "{}".format(Lam_w)\
+"_lambda_nw" + "{}".format(Lam_nw)\
+"_Lw" + "{}".format(L["wetting"])\
- +"_Lnw" + "{}".format(L["nonwetting"])\
- +".xdmf"
+ +"_Lnw" + "{}".format(L["nonwetting"])
+ filename = self.output_dir+"_solution" +".xdmf"
# create solution files and populate dictionary.
+ self.output_filename_parameter_part.update(
+ {subdom_ind: filename_param_part}
+ )
self.solution_file.update( #
{subdom_ind: SolutionFile(self._mpi_communicator, filename)}
)
- def write_subsequent_errors(self,#
+ def write_subsequent_errors_to_csv(self,#
filename: str, #
subdomain_index: int,#
- subsequent_errors,#: tp.Dict[str: float],#
+ errors: tp.Dict[str, float],#
time: float = None) -> None:
""" write subsequent errors to csv file for plotting with pgfplots"""
# in case no time is given, use the time of the class.
@@ -565,14 +615,48 @@ class LDDsimulation(object):
self.errors = pd.DataFrame({ "iteration": iteration,
#"th. initial mass": df.assemble(rho*df.dx), \
#"energy": self.total_energy(), \
- "wetting": subsequent_errors["wetting"]}, index=[iteration])
+ "wetting": errors["wetting"]}, index=[iteration])
else:
self.errors = pd.DataFrame({"iteration":iteration, #\
#"th. initial mass": df.assemble(rho*df.dx), \
#"energy": self.total_energy(), \
- "wetting": subsequent_errors["wetting"],
- "nonwetting": subsequent_errors["nonwetting"]}, index=[iteration])
+ "wetting": errors["wetting"],
+ "nonwetting": errors["nonwetting"]}, index=[iteration])
+ # if self._rank == 0:
+ # check if file exists
+ if os.path.isfile(filename) == True:
+ with open(filename, 'a') as f:
+ self.errors.to_csv(f, header=False, sep='\t', encoding='utf-8', index=False)
+ else:
+ self.errors.to_csv(filename, sep='\t', encoding='utf-8', index=False)
+
+ def write_errornorms_to_csv(self,#
+ filename: str, #
+ subdomain_index: int,#
+ errors: tp.Dict[str, float],#
+ ) -> 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,
+ #"th. initial mass": df.assemble(rho*df.dx), \
+ #"energy": self.total_energy(), \
+ "wetting": errors["wetting"]}, index=[timestep])
+
+ else:
+ self.errors = pd.DataFrame({ "time": self.t,
+ "timestep": self.timestep_num,
+ "wetting": errors["wetting"],
+ "nonwetting": errors["nonwetting"]}, index=[timestep])
# if self._rank == 0:
# check if file exists
if os.path.isfile(filename) == True:
@@ -866,7 +950,7 @@ class LDDsimulation(object):
# subdomain.outer_boundary[j]['nonwetting'] return
# the actual expression needed for the dirichlet condition for both phases if present.
for index, boundary in enumerate(subdomain.outer_boundary):
- if np.abs(time) < self.tol:
+ if np.abs(time-self.t_initial) < self.tol:
# Here the dictionary has to be created first because t = 0.
# create a dictionary of dictionaries for each subdomain
self.outerBC.update({num: dict()})
@@ -885,7 +969,7 @@ class LDDsimulation(object):
boundary_marker_value = index+1
for phase in subdomain.has_phases:
# note that the default interpolation degree is 2
- if np.abs(time) < self.tol :
+ if np.abs(time - self.t_initial) < self.tol :
# time = 0 so the Dolfin Expression has to be created first.
pDCa = df.Expression(pDC[phase], domain = mesh, degree = interpolation_degree, t=time)
self.dirichletBC_dfExpression[num][index].update({phase: pDCa})
@@ -900,6 +984,27 @@ class LDDsimulation(object):
print("subdomain is an inner subdomain and has no outer boundary.\n",
"no Dirichlet boundary has been set.")
+ def _init_exact_solution_expression(self):
+ """ initialise dolfin expressions for exact solution and populate
+ self.exact_solution_expr dictionary. This is used to calculate errornorms
+ in the case of an exact solution.
+ """
+ if self.exact_solution is not None:
+ for sd_ind, subdomain in self.subdomain.items():
+ for phase in subdomain.has_phases:
+ if np.abs(self.t-self.t_initial) < self.calc_tol:
+ # here t = 0 and we have to initialise the sources.
+ mesh = subdomain.mesh
+ # p0 contains both pw_0 and pnw_0 for subdomain[subdomain_index]
+ pa_exact = self.exact_solution[sd_ind][phase]
+ self.exact_solution_expr.update({sd_ind: dict()})
+ # note that the default interpolation degree is 2
+ self.exact_solution_expr[sd_ind].update(
+ { phase: df.Expression(pa_exact, domain = mesh, degree=self.interpolation_degree, t=self.t) }
+ )
+ else:
+ pass
+
def _eval_sources(self, subdomain_index: int, #
interpolation_degree: int = None, #
time: float = 0,#
@@ -911,7 +1016,7 @@ class LDDsimulation(object):
subdomain = self.subdomain[subdomain_index]
for phase in subdomain.has_phases:
- if np.abs(time) < self.calc_tol:
+ if np.abs(time-self.t_initial) < self.calc_tol:
# here t = 0 and we have to initialise the sources.
mesh = subdomain.mesh
# p0 contains both pw_0 and pnw_0 for subdomain[subdomain_index]
@@ -928,6 +1033,18 @@ class LDDsimulation(object):
# note that the default interpolation degree is 2
subdomain.source[phase].t = time
+ def _add_parameters_to_output_dirname(self):
+ """ function to create a new output directory for given set of input
+ parameters and adjust self.output_dir accordingly
+ """
+ dirname = "mesh_res" + "{}_".format(self.mesh_resolution)\
+ + "dt" + "{}".format(self.timestep_size) + "/"
+ self.output_dir += dirname
+
+ if not os.path.exists(self.output_dir):
+ os.mkdir(self.output_dir)
+ print(f"\nwill save output of this simulation to path {self.output_dir}")
+
# https://github.com/geo-fluid-dynamics/phaseflow-fenics
# adapted from Lucas Ostrowski
diff --git a/LDDsimulation/domainPatch.py b/LDDsimulation/domainPatch.py
index deae328..8b3ce51 100644
--- a/LDDsimulation/domainPatch.py
+++ b/LDDsimulation/domainPatch.py
@@ -288,7 +288,7 @@ class DomainPatch(df.SubDomain):
interface_forms = []
# the marker values of the interfacese are global interfaces index + 1
for interface in self.has_interface:
- print(f"subdomain{self.subdomain_index} has interfaces: {interface}\n")
+ # print(f"subdomain{self.subdomain_index} has interfaces: {interface}\n")
marker_value = self.interface[interface].marker_value
# interface_forms += (dt*Lambda_a*pa*phi_a)*ds(interface)
interface_forms.append((dt*Lambda_a*pa*phi_a)*ds(marker_value))
@@ -665,18 +665,18 @@ class DomainPatch(df.SubDomain):
mesh_data = self.mesh.data()
# local to parent vertex index map
self.parent_mesh_index = mesh_data.array('parent_vertex_indices',0)
- print(f"on subdomain{self.subdomain_index} parent_vertex_indices: {self.parent_mesh_index}")
+ # print(f"on subdomain{self.subdomain_index} parent_vertex_indices: {self.parent_mesh_index}")
self.dof2vertex = dict()
self.vertex2dof = dict()
for phase in self.has_phases:
self.dof2vertex.update(#
{phase :df.dof_to_vertex_map(self.function_space[phase])}#
)
- print(f"on subdomain{self.subdomain_index} dof2vertex: {self.dof2vertex[phase]}")
+ # print(f"on subdomain{self.subdomain_index} dof2vertex: {self.dof2vertex[phase]}")
self.vertex2dof.update(#
{phase :df.vertex_to_dof_map(self.function_space[phase])}#
)
- print(f"on subdomain{self.subdomain_index} vertex2dof: {self.vertex2dof[phase]}")
+ # print(f"on subdomain{self.subdomain_index} vertex2dof: {self.vertex2dof[phase]}")
def _init_markers(self) -> None:
""" define boundary markers
@@ -1124,13 +1124,13 @@ class Interface(BoundaryPart):
"""
interface_vertex_indices = self._vertex_indices(interface_marker, interface_marker_value,
print_vertex_indices = False)
- print(f"local interface vertices: {interface_vertex_indices}")
+ # print(f"local interface vertices: {interface_vertex_indices}")
dof2vertex = df.dof_to_vertex_map(function_space)
- print(f"dof2vertex map: {dof2vertex}\n")
+ # print(f"dof2vertex map: {dof2vertex}\n")
# check which dofs actually lie on the interface using the dof2vertex map.
is_a_dof_on_interface = np.isin(dof2vertex, interface_vertex_indices, assume_unique=True)
- print(f"is_a_dof_on_interface: {is_a_dof_on_interface}")
- print(f"this means the following dofs are interface dofs:", np.nonzero(is_a_dof_on_interface)[0])
+ # print(f"is_a_dof_on_interface: {is_a_dof_on_interface}")
+ # print(f"this means the following dofs are interface dofs:", np.nonzero(is_a_dof_on_interface)[0])
return np.nonzero(is_a_dof_on_interface)[0]
def init_interface_values(self,#
diff --git a/RR-2-patch-test-case/RR-2-patch-test.py b/RR-2-patch-test-case/RR-2-patch-test.py
index c61dff9..5857fe0 100755
--- a/RR-2-patch-test-case/RR-2-patch-test.py
+++ b/RR-2-patch-test-case/RR-2-patch-test.py
@@ -84,10 +84,40 @@ isRichards = {
2: True
}
-Tmax = 1
-timestep_size = 0.005
+##################### MESH #####################################################
+mesh_resolution = 40
+##################### TIME #####################################################
+timestep_size = 2*0.0001
+number_of_timesteps = 300
+# 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_analise = 7
+starttime = 0
+
+# We want to write out csv files containing subsequent L2 errors of iterations for
+# number_of_timesteps_to_analise many timesteps. We want to do this with
+# analyse_timesteps = np.linspace(starttime, number_of_timesteps, number_of_timesteps_to_analise)
+# see below. Therefor, we need to check if (number_of_timesteps_to_analise-1) is
+# a divisor of number_of_timesteps,
+if number_of_timesteps_to_analise == 0:
+ analyse_timesteps = None
+elif number_of_timesteps_to_analise == 1:
+ analyse_timesteps = [0]
+else:
+ if not number_of_timesteps%(number_of_timesteps_to_analise-1) == 0:
+ # if number_of_timesteps%(number_of_timesteps_to_analise-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_analise = {number_of_timesteps/number_of_timesteps_to_analise}",
+ "is not an integer,\n")
+ new_analysing_interval = int(round(number_of_timesteps/(number_of_timesteps_to_analise-1)))
+ number_of_timesteps = new_analysing_interval*(number_of_timesteps_to_analise-1)
+ print(f"I'm resetting number_of_timesteps = {number_of_timesteps}\n")
+ analyse_timesteps = np.linspace(0, number_of_timesteps, number_of_timesteps_to_analise)
+
+Tmax = number_of_timesteps*timestep_size
+time_interval = [starttime, Tmax]
+print(f"The simulation interval is [{starttime}, {Tmax}]")
-time_interval = [0, Tmax]
viscosity = {#
# subdom_num : viscosity
@@ -109,8 +139,8 @@ L = {#
lambda_param = {#
# subdom_num : lambda parameter for the L-scheme
- 1 : {'wetting' :100},#
- 2 : {'wetting' :100}
+ 1 : {'wetting' :160},#
+ 2 : {'wetting' :160}
}
## relative permeabilty functions on subdomain 1
@@ -190,8 +220,6 @@ write_to_file = {
'L_iterations': True
}
-mesh_resolution = 30
-
# initialise LDD simulation class
simulation = ldd.LDDsimulation(tol = 1E-12)
simulation.set_parameters(output_dir = "./output/",#
@@ -212,348 +240,11 @@ simulation.set_parameters(output_dir = "./output/",#
sources = source_expression,#
initial_conditions = initial_condition,#
dirichletBC_Expressions = dirichletBC,#
- exact_solution_case = True,#
+ exact_solution = exact_solution,#
write2file = write_to_file,#
)
simulation.initialise()
-# simulation._init_meshes_and_markers(subdomain_def_points, mesh_resolution=2)
-# subdomain marker functions
-# output_dir = simulation.output_dir
-# domain_marker = simulation.domain_marker
-# mesh_subdomain = simulation.mesh_subdomain
-
-
-# mesh = mesh_subdomain[0]
-# interface_marker = df.MeshFunction('int', mesh, mesh.topology().dim()-1)
-# interface_marker.set_all(0)
-# interface = dp.Interface(vertices=interface_def_points[0], #
-# tol=mesh.hmin()/100,#
-# internal=True,#
-# adjacent_subdomains = adjacent_subdomains[0])
-# print("interface vertices are: \n", interface.vertices)
-# interface.mark(interface_marker, 1)
-# simulation._init_interfaces(interface_def_points, adjacent_subdomains)
-
-# interface = simulation.interface
-# interface_marker = simulation.interface_marker
-#
-# subdoms = simulation.subdomain
-
-# df.File('./subdomain1_interface_marker.pvd') << simulation.subdomain[1].interface_marker
-# df.File('./subdomain2_interface_marker.pvd') << simulation.subdomain[2].interface_marker
-# df.File('./test_subdomain1_boundary_marker.pvd') << boundary_marker1
-# df.File('./test_subdomain2_boundary_marker.pvd') << boundary_marker2
-analyse_timesteps = [0, 100*timestep_size, 1, 2, 3]
simulation.run(analyse_timesteps)
# simulation.LDDsolver(time = 0, debug = True, analyse_timestep = True)
# df.info(parameters, True)
-
-# for iterations in range(1):
-# #
-#
-# a1 = (L1*u1*v1)*dx1 + (df.inner(relative_permeability(saturation(p1_i, 1), 1)*df.grad(u1), df.grad(v1)))*dx_1 + (timestep*lambda1*u1*v1)*ds1(1) #timestep*
-# rhs1 = (L1*p1_i*v1)*dx1 - ((saturation(p1_i, 1) - saturation(p1_0, 1))*v1)*dx1 + (timestep*(source1 - g1_i)*v1)*ds1(1)
-# rhssplit1 = (L1*p1_i*v1)*dx1 - ((saturation(p1_i, 1) - saturation(p1_0, 1))*v1)*dx1
-# extratrem = (timestep*(source1 - g1_i)*v1)*ds1(1)
-# splitrhs = df.assemble(rhssplit1)
-# print("splitrhs: \n", splitrhs.get_local())
-# extrat_assembled = df.assemble(extratrem)
-# print("extratrem: \n", extrat_assembled.get_local())
-# added = splitrhs + extrat_assembled
-# print("both added together:\n", added.get_local())
-#
-# A1 = df.assemble(a1)
-# dsform = (timestep*lambda1*u1*v1)*ds1(1)
-# dsform_vec = df.assemble(dsform)
-# b1 = df.assemble(rhs1)
-# #p_gamma1.apply(A1,b1)
-# # outerBC1.apply(A1,b1)
-# u = df.Function(V1)
-# U1 = u.vector()
-# df.solve(A1,U1,b1)
-# # print("Form:\n", A1.array())
-# print("RHS:\n", b1.get_local())
-# #print("ds term:\n", dsform_vec.array())
-# # print('solution:\n', U1.get_local())
-# p1_i.vector()[:] = U1
-# interface[0].read_pressure_dofs(from_function = p1_i, #
-# interface_dofs = dofs_on_interface1,#
-# dof_to_vert_map = dom1_d2v,#
-# local_to_parent_vertex_map = dom1_parent_vertex_indices,#
-# phase = 'wetting',#
-# subdomain_ind = 1)
-# #
-# # Save mesh to file
-# df.File('./test_domain_layered_soil.xml.gz') << mesh_subdomain[0]
-# df.File('./test_domain_mesh.pvd') << mesh_subdomain[0]
-# df.File('./test_global_interface_marker.pvd') << interface_marker
-# df.File('./test_subdomain1.xml.gz') << mesh_subdomain[1]
-# df.File('./test_subdomain2.xml.gz') << mesh_subdomain[2]
-# df.File('./test_domain_marker.pvd') << domain_marker
-# df.File('./test_domain_layered_soil_solution.pvd') << u
-# df.File('./test_subdomain1_interface_marker.pvd') << interface_marker1
-# df.File('./test_subdomain2_interface_marker.pvd') << interface_marker2
-# df.File('./test_subdomain1_boundary_marker.pvd') << boundary_marker1
-# df.File('./test_subdomain2_boundary_marker.pvd') << boundary_marker2
-
-# boundary_marker1 = simulation.subdomain[1].outer_boundary_marker
-# boundary_marker2 = simulation.subdomain[2].outer_boundary_marker
-# dx_1 = simulation.subdomain[1].dx
-# dx_2 = simulation.subdomain[2].dx
-# ds_1 = simulation.subdomain[1].ds
-# ds_2 = simulation.subdomain[2].ds
-#
-# interface_marker1 = simulation.subdomain[1].interface_marker#df.MeshFunction('int', mesh_subdomain[1], mesh_subdomain[1].topology().dim()-1)
-# interface_marker2 = simulation.subdomain[2].interface_marker#df.MeshFunction('int', mesh_subdomain[2], mesh_subdomain[2].topology().dim()-1)
-# # interface_marker1.set_all(0)
-# # interface_marker2.set_all(0)
-# # print(dir(interface_marker1))
-#
-# interface[0].mark(interface_marker1, 1)
-# interface[0].mark(interface_marker2, 1)
-# # interface_marker1 = simulation.subdomain[1].boundary_marker
-#
-# # domain measures
-#
-# # dx1 = simulation.subdomain[1].dx
-# dx1 = df.Measure('dx', domain = mesh_subdomain[1])
-# dx2 = df.Measure('dx', domain = mesh_subdomain[2])
-#
-# # boundary measures
-# ds1 = df.Measure('ds', domain = mesh_subdomain[1], subdomain_data = interface_marker1)
-# ds2 = df.Measure('ds', domain = mesh_subdomain[2], subdomain_data = interface_marker2)
-#
-# V1 = df.FunctionSpace(mesh_subdomain[1], 'P', 1)
-# p1_exact = df.Expression('1.0 - (1.0 + t*t)*(1.0 + x[0]*x[0] + x[1]*x[1])', domain = mesh_subdomain[1], degree = 2, t = 0)
-#
-# coordinates = mesh_subdomain[0].coordinates()
-# print("global mesh coordinates of domain: \n", coordinates)
-# coordinates1 = mesh_subdomain[1].coordinates()
-#
-# coordinates2 = mesh_subdomain[2].coordinates()
-# submesh1_data = mesh_subdomain[1].data()
-# dom1_parent_vertex_indices = submesh1_data.array('parent_vertex_indices',0)
-# # print("map of dom1 vertices to parent vertices on global domain:\n", dom1_parent_vertex_indices)
-# #
-# # coordinates2 = mesh_subdomain[2].coordinates()
-# # print("\nvertices of subdomain2: \n", coordinates2)
-# submesh2_data = mesh_subdomain[2].data()
-# dom2_parent_vertex_indices = submesh2_data.array('parent_vertex_indices',0)
-# # print("'parent_vertex_indices',1 \n")
-# #
-# #
-# # print("on domain: \n")
-# # interface_coordinates = interface[0].coordinates(interface_marker, 1, print_coordinates = True)
-# # interface_coordinates = interface[0]._vertex_indices(interface_marker, 1, print_vertex_indices = True)
-# # print("\non subdomain2: \n")
-# # dom2_interface_coordinates = interface[0].coordinates(interface_marker2, 1, print_coordinates = True)
-# dom2_interface_def_points = interface[0]._vertex_indices(interface_marker2, 1, print_vertex_indices = True)
-#
-#
-#
-# dofs_on_interface1 = interface[0].dofs_on_interface(V1, interface_marker1, 1)
-# print("Dofs on Interface dom1", dofs_on_interface1)
-#
-# V2 = df.FunctionSpace(mesh_subdomain[2], 'P', 1)
-# dofs_on_interface2 = interface[0].dofs_on_interface(V2, interface_marker2, 1)
-# # print("Dofs on Interface dom1", dofs_on_interface2)
-#
-# # markermesh = interface_marker2.mesh()
-# # functionSpaceMesh = V2.mesh()
-# # print("Domain 2 mesh extracted from markerfunction:\n", markermesh.coordinates())
-# # print("Domain 2 mesh extracted from Function Space V2:\n", functionSpaceMesh.coordinates())
-#
-# dom1_d2v = df.dof_to_vertex_map(V1)
-# dom1_v2d = df.vertex_to_dof_map(V1)
-# # print("dof to vertex map for V1:\n", dom1_d2v)
-# # print("vertex to dof map for V1:\n", dom1_v2d)
-# dom2_d2v = df.dof_to_vertex_map(V2)
-# dom2_v2d = df.vertex_to_dof_map(V2)
-# # print("dof to vertex map for V2:\n", dom2_d2v)
-# # print("vertex to dof map for V2:\n", dom2_v2d)
-#
-# testf1 = df.Constant(41)
-# f1test = df.interpolate(testf1, V1)
-# testf2 = df.Constant(42)
-# f2test = df.interpolate(testf2, V2)
-#
-# # f1vec = f1test.compute_vertex_values()
-# # for i in range(len(f1vec)):
-# # f1vec[i] = i
-#
-# # f1test.vector()[:] = f1vec
-# # f1nodal = f1test.vector().get_local()
-# # for i in range(len(f1nodal)):
-# # f1nodal[i] = i
-#
-# # print("nodal values of f1nodal\n", f1nodal)
-#
-# #
-# # print('dofs on interface \n', f1vec[dofs_on_interface])
-# # f1test.vector()[:] = f1vec
-# # print("f1test.vector()[:]\n", f1test.vector()[:])
-# # print("f1test.vector()[dofs_on_interface]\n", f1test.vector()[dofs_on_interface])
-#
-# # print("constant function on dom1\n", f1vec)
-# #
-# # for i in dofs_on_interface1:
-# # print("local dof ind %i, parent = %i", i,dom1_parent_vertex_indices[dom1_d2v[i]])
-# # # f1test.vector() is enumerated by dof not by vertex
-# # f1test.vector()[i] = i
-# # # this enumerates according to the dof numbering
-# # print("testfunction f1test with set interface values\n", f1test.vector()[:])
-# # print("testfunction f1test with get_local\n", f1test.vector().get_local())
-# # f1test.vector()[dom1_v2d[dom1_d2v[dofs_on_interface]]] = f1vec[dofs_on_interface]
-# # .vector()[:] and .vector().get_local() order values according to the vertex numbering
-# # print("testfunction f1test with compute_vertex_values\n", f1test.compute_vertex_values())
-# # print("testfunction f1test with compute_vertex_values mit d2vmap\n", f1test.compute_vertex_values()[dom1_d2v[:]])
-# # i=0
-# # for x in coordinates1:
-# # print(f"vertex {i}: f1test[{dom1_v2d[i]}] = {f1test.vector()[dom1_v2d[i]]}\tu({x}) = {f1test(x)}")
-# # i += 1
-# # print("\non subdomain1: \n")
-# print("vertices of subdomain1: \n", coordinates1)
-# # dom1_interface_coordinates = interface[0].coordinates(interface_marker1, 1, print_coordinates = True)
-# dom1_interface_def_points = interface[0]._vertex_indices(interface_marker1, 1, print_vertex_indices = True)
-#
-# interface[0].read_pressure_dofs(from_function = f1test, #
-# interface_dofs = dofs_on_interface1,#
-# dof_to_vert_map = dom1_d2v,#
-# local_to_parent_vertex_map = dom1_parent_vertex_indices,#
-# phase = 'wetting',
-# subdomain_ind = 1)
-#
-# interface[0].write_pressure_dofs(to_function = f2test, #
-# interface_dofs = dofs_on_interface2,#
-# dof_to_vert_map = dom2_d2v,#
-# local_to_parent_vertex_map = dom2_parent_vertex_indices,#
-# phase = 'wetting',
-# subdomain_ind = 2)
-#
-# i=0
-# for x in coordinates2:
-# print(f"vertex {i}: f2test[{dom2_v2d[i]}] = {f2test.vector()[dom2_v2d[i]]}\tu({x}) = {f2test(x)}")
-# i += 1
-# print("\non subdomain2: \n")
-# print("vertices of subdomain2: \n", coordinates2)
-#
-# # dom2_interface_coordinates = interface[0].coordinates(interface_marker2, 1, print_coordinates = True)
-# dom2_interface_def_points = interface[0]._vertex_indices(interface_marker2, 1, print_vertex_indices = True)
-# print("Parent indices", dom2_parent_vertex_indices)
-#
-#
-# u1 = df.TrialFunction(V1)
-# v1 = df.TestFunction(V1)
-#
-# u2 = df.TrialFunction(V2)
-# v2 = df.TestFunction(V2)
-#
-# def relative_permeability(s, subdomain_index):
-# if subdomain_index == 1:
-# # relative permeabilty on subdomain1
-# return s**2
-# else:
-# # relative permeabilty on subdomain2
-# return s**3
-#
-# def saturation(pressure, subdomain_index):
-# # inverse capillary pressure-saturation-relationship
-# return df.conditional(pressure < 0, 1/((1 - pressure)**(1/(subdomain_index + 1))), 1)
-#
-#
-# # exact solution
-# # cell = df.Cell("triangle", 2)
-# # x1 = df.SpatialCoordinate(mesh_subdomain[1])
-# # x2 = df.SpatialCoordinate(mesh_subdomain[2])
-#
-# p2_exact = df.Expression('1.0 - (1.0 + t*t)*(1.0 + x[1]*x[1])', domain = mesh_subdomain[2], degree = 2, t = 0)
-# p1_e = df.interpolate(p1_exact, V1)
-# p2_e = df.interpolate(p2_exact, V2)
-# #
-# # initial conditions
-# p1_initial = df.Expression('-(x[0]*x[0] + x[1]*x[1])', domain = mesh_subdomain[1], degree = 2)
-# p2_initial = df.Expression('-x[1]*x[1]', domain = mesh_subdomain[2], degree = 2)
-# p1_0 = df.interpolate(p1_initial, V1)
-# p2_0 = df.interpolate(p2_initial, V2)
-#
-# # Initial boundary value for the pressure on interface
-# p_gamma1 = df.DirichletBC(V1, p1_exact, interface_marker1, 1)
-# p_gamma2 = df.DirichletBC(V2, p2_exact, interface_marker2, 1)
-#
-# outerBC1 = df.DirichletBC(V1, p1_exact, boundary_marker1, 1)
-# outerBC2 = df.DirichletBC(V2, p2_exact, boundary_marker2, 1)
-#
-#
-# ### source terms
-# # source1 = df.Expression('4.0/pow(1 + x[0]*x[0] + x[1]*x[1], 2) - t/sqrt( pow(1 + t*t, 3)*(1 + x[0]*x[0] + x[1]*x[1]) )', domain = mesh_subdomain[1], degree = 2, t = 0)
-# # source2 = df.Expression('2.0*(1-x[1]*x[1])/pow(1 + x[1]*x[1], 2) - 2*t/(3*pow( pow(1 + t*t, 4)*(1 + x[1]*x[1]), 1/3))', domain = mesh_subdomain[2], degree = 2, t = 0)
-# source1 = df.Expression(source_expression[1]['wetting'], domain = mesh_subdomain[1], degree = 2, t = 0)
-# source2 = df.Expression('2.0*(1-x[1]*x[1])/pow(1 + x[1]*x[1], 2) - 2*t/(3*pow( pow(1 + t*t, 4)*(1 + x[1]*x[1]), 1/3))', domain = mesh_subdomain[2], degree = 2, t = 0)
-#
-# source1h = df.interpolate(source1, V1)
-# source2h = df.interpolate(source2, V2)
-#
-# L1 = 0.25
-# L2 = L1
-# timestep = 0.1
-# lambda1 = 4
-# lambda2 = lambda1
-# # initial iteartion
-# p1_i = p1_0
-# p2_i = p2_0
-# # initial g_i
-# n1 = df.FacetNormal(mesh_subdomain[1])
-# flux = -df.grad(p1_0)
-# gl0 = (df.dot(flux, n1) - lambda1*p1_0)*v1*ds1(1)
-# print('gl0 term assembliert', df.assemble(gl0).get_local())
-# g1_i = -lambda1*p1_i
-# g2_i = -lambda1*p2_i
-
-# for iterations in range(1):
-# #
-#
-# a1 = (L1*u1*v1)*dx1 + (df.inner(relative_permeability(saturation(p1_i, 1), 1)*df.grad(u1), df.grad(v1)))*dx_1 + (timestep*lambda1*u1*v1)*ds1(1) #timestep*
-# rhs1 = (L1*p1_i*v1)*dx1 - ((saturation(p1_i, 1) - saturation(p1_0, 1))*v1)*dx1 + (timestep*(source1 - g1_i)*v1)*ds1(1)
-# rhssplit1 = (L1*p1_i*v1)*dx1 - ((saturation(p1_i, 1) - saturation(p1_0, 1))*v1)*dx1
-# extratrem = (timestep*(source1 - g1_i)*v1)*ds1(1)
-# splitrhs = df.assemble(rhssplit1)
-# print("splitrhs: \n", splitrhs.get_local())
-# extrat_assembled = df.assemble(extratrem)
-# print("extratrem: \n", extrat_assembled.get_local())
-# added = splitrhs + extrat_assembled
-# print("both added together:\n", added.get_local())
-#
-# A1 = df.assemble(a1)
-# dsform = (timestep*lambda1*u1*v1)*ds1(1)
-# dsform_vec = df.assemble(dsform)
-# b1 = df.assemble(rhs1)
-# #p_gamma1.apply(A1,b1)
-# # outerBC1.apply(A1,b1)
-# u = df.Function(V1)
-# U1 = u.vector()
-# df.solve(A1,U1,b1)
-# # print("Form:\n", A1.array())
-# print("RHS:\n", b1.get_local())
-# #print("ds term:\n", dsform_vec.array())
-# # print('solution:\n', U1.get_local())
-# p1_i.vector()[:] = U1
-# interface[0].read_pressure_dofs(from_function = p1_i, #
-# interface_dofs = dofs_on_interface1,#
-# dof_to_vert_map = dom1_d2v,#
-# local_to_parent_vertex_map = dom1_parent_vertex_indices,#
-# phase = 'wetting',#
-# subdomain_ind = 1)
-# #
-# # Save mesh to file
-# df.File('./test_domain_layered_soil.xml.gz') << mesh_subdomain[0]
-# df.File('./test_domain_mesh.pvd') << mesh_subdomain[0]
-# df.File('./test_global_interface_marker.pvd') << interface_marker
-# df.File('./test_subdomain1.xml.gz') << mesh_subdomain[1]
-# df.File('./test_subdomain2.xml.gz') << mesh_subdomain[2]
-# df.File('./test_domain_marker.pvd') << domain_marker
-# df.File('./test_domain_layered_soil_solution.pvd') << u
-# df.File('./test_subdomain1_interface_marker.pvd') << interface_marker1
-# df.File('./test_subdomain2_interface_marker.pvd') << interface_marker2
-# df.File('./test_subdomain1_boundary_marker.pvd') << boundary_marker1
-# df.File('./test_subdomain2_boundary_marker.pvd') << boundary_marker2
diff --git a/RR-2-patch-test-case/plots/exact_solution_error_norms.tex b/RR-2-patch-test-case/plots/exact_solution_error_norms.tex
new file mode 100644
index 0000000..8d7d343
--- /dev/null
+++ b/RR-2-patch-test-case/plots/exact_solution_error_norms.tex
@@ -0,0 +1,47 @@
+\documentclass[a4paper]{scrartcl}%scrartcl
+\input{praeambels_and_definitions/artikelpraeambel.tex}
+\input{praeambels_and_definitions/artikel_theoreme_und_farbe.tex}
+%\input{praeambels_and_definitions/special_symbols.tex}
+%%%% MAKROS %%%%%%%%%%%%% restriction
+\input{./praeambels_and_definitions/L-Schema_Paper_makros.tex}
+% set graphicspath
+\graphicspath{ {./} }
+
+\usepackage{pgfplotstable}
+
+\begin{document}
+ % Generation of Profile_Error_new.png
+ \tikzsetnextfilename{Errornorm-exact-solution-over-time}
+ \tikzset{external/force remake}
+ \begin{tikzpicture}
+% % \path[draw,dashed,thick] (3.13,0) -- (3.13,5.2);
+ \begin{axis}[%
+ width=\textwidth,
+ title={ $\bigl\|p_{w,l} - p_{w,l}^n\bigr\|_{L^2}$ over time $t$, $ h \approx 0.02$, $\tau = 2\cdot 10^{-4}$ },
+ % axis lines=left,
+ % legend style = {draw=none},
+ legend cell align = left,
+ xlabel= {timesteps},
+ ylabel= {$\bigl\|p_{w,l} - p_{w,l}^n\bigr\|_{L^2}$},
+% xmin= 0,
+% xmax= 53,
+% ymin= 0,
+% ymax= 0.0003,
+ grid= both, %major or minor
+ axis line style={-Latex[round]},
+ legend style={
+% anchor=north east,
+% at={(1,1)},
+ font=\tiny
+ },
+ legend entries={$\bigl\|p_{w,1} - p_{w,1}^{n}\bigr\|_{L^2(\dom_1)}$,
+ $\bigl\|p_{w,2} - p_{w,2}^{n}\bigr\|_{L^2(\dom_2)}$},
+ ] % end of axis options %col sep=comma,
+ \addplot table [x=timestep, y=wetting] {../output/mesh_res30_dt0.004/subdomain1_dt0.004_hmin0.019_lambda70_Lw0.25_L2_errornorms_over_time.csv};%
+ \addplot table [x=timestep, y=wetting] {../output/mesh_res30_dt0.004/subdomain2_dt0.004_hmin0.019_lambda70_Lw0.25_L2_errornorms_over_time.csv};%
+% % \addplot table [col sep=comma] {../output/subdomain1_subsequent_iteration_error_for_phase_wetting_at_time0.csv};%
+% % % \addplot table [col sep=comma] {../output/subdomain2_subsequent_iteration_error_for_phase_wetting_at_time0.csv};%
+ \end{axis}
+ \end{tikzpicture}
+
+\end{document}
diff --git a/RR-2-patch-test-case/plots/subsequent_errors.tex b/RR-2-patch-test-case/plots/subsequent_errors.tex
index 2d98427..7df36b5 100644
--- a/RR-2-patch-test-case/plots/subsequent_errors.tex
+++ b/RR-2-patch-test-case/plots/subsequent_errors.tex
@@ -10,21 +10,21 @@
\usepackage{pgfplotstable}
\begin{document}
- % Generation of Profile_Error_new.png
- \tikzsetnextfilename{Subsequent_errors_t0}
+ % Generation of Subsequent_errors_t0.00
+ \tikzsetnextfilename{Subsequent_errors_t0.00}
\tikzset{external/force remake}
\begin{tikzpicture}
% % \path[draw,dashed,thick] (3.13,0) -- (3.13,5.2);
\begin{semilogyaxis}[%
width=\textwidth,
- title={ Subsequent errors for $t = 0$, $ h \approx 0.02$, $\tau = 2\cdot 10^{-4}$ },
+ title={ Subsequent errors for $t = 0$, $ h \approx 0.02$, $\tau = 4\cdot 10^{-4}$ },
% axis lines=left,
% legend style = {draw=none},
legend cell align = left,
xlabel= {iterations},
ylabel= {subsequent errors},
xmin= 0,
- xmax= 53,
+% xmax= 53,
% ymin= 0,
% ymax= 0.0003,
grid= both, %major or minor
@@ -37,8 +37,110 @@
legend entries={$\bigl\|p_l^i - p_l^{i-1}\bigr\|_{L^2(\dom_1)}$,
$\bigl\|p_l^i - p_l^{i-1}\bigr\|_{L^2(\dom_2)}$},
] % end of axis options %col sep=comma,
- \addplot table [x=iteration, y=wetting] {../output/subdomain1_subsequent_iteration_error_for_phase_wetting_at_time0.csv};%
- \addplot table [x=iteration, y=wetting] {../output/subdomain2_subsequent_iteration_error_for_phase_wetting_at_time0.csv};%
+ \addplot table [x=iteration, y=wetting] {../output/mesh_res30_dt0.004/subdomain1_dt0.004_hmin0.019_lambda70_Lw0.25subsequent_iteration_errors_at_time0.00.csv};%
+ \addplot table [x=iteration, y=wetting] {../output/mesh_res30_dt0.004/subdomain2_dt0.004_hmin0.019_lambda70_Lw0.25subsequent_iteration_errors_at_time0.00.csv};%
+% \addplot table [col sep=comma] {../output/subdomain1_subsequent_iteration_error_for_phase_wetting_at_time0.csv};%
+% % \addplot table [col sep=comma] {../output/subdomain2_subsequent_iteration_error_for_phase_wetting_at_time0.csv};%
+ \end{semilogyaxis}
+ \end{tikzpicture}
+
+ % Generation of Subsequent_errors_t0.00
+ \tikzsetnextfilename{Subsequent_errors_t0.05}
+ \tikzset{external/force remake}
+ \begin{tikzpicture}
+% % \path[draw,dashed,thick] (3.13,0) -- (3.13,5.2);
+ \begin{semilogyaxis}[%
+ width=\textwidth,
+ title={ Subsequent errors for $t = 0.05$, $ h \approx 0.02$, $\tau = 4\cdot 10^{-4}$ },
+ % axis lines=left,
+ % legend style = {draw=none},
+ legend cell align = left,
+ xlabel= {iterations},
+ ylabel= {subsequent errors},
+% xmin= 0,
+% xmax= 53,
+% ymin= 0,
+% ymax= 0.0003,
+ grid= both, %major or minor
+ axis line style={-Latex[round]},
+ legend style={
+% anchor=north east,
+% at={(1,1)},
+ font=\tiny
+ },
+ legend entries={$\bigl\|p_l^i - p_l^{i-1}\bigr\|_{L^2(\dom_1)}$,
+ $\bigl\|p_l^i - p_l^{i-1}\bigr\|_{L^2(\dom_2)}$},
+ ] % end of axis options %col sep=comma,
+ \addplot table [x=iteration, y=wetting] {../output/mesh_res30_dt0.004/subdomain1_dt0.004_hmin0.019_lambda70_Lw0.25subsequent_iteration_errors_at_time0.05.csv};%
+ \addplot table [x=iteration, y=wetting] {../output/mesh_res30_dt0.004/subdomain2_dt0.004_hmin0.019_lambda70_Lw0.25subsequent_iteration_errors_at_time0.05.csv};%
+% \addplot table [col sep=comma] {../output/subdomain1_subsequent_iteration_error_for_phase_wetting_at_time0.csv};%
+% % \addplot table [col sep=comma] {../output/subdomain2_subsequent_iteration_error_for_phase_wetting_at_time0.csv};%
+ \end{semilogyaxis}
+ \end{tikzpicture}
+
+ % Generation of Subsequent_errors_t0.00
+ \tikzsetnextfilename{Subsequent_errors_t0.1}
+ \tikzset{external/force remake}
+ \begin{tikzpicture}
+% % \path[draw,dashed,thick] (3.13,0) -- (3.13,5.2);
+ \begin{semilogyaxis}[%
+ width=\textwidth,
+ title={ Subsequent errors for $t = 0.1$, $ h \approx 0.02$, $\tau = 4\cdot 10^{-4}$ },
+ % axis lines=left,
+ % legend style = {draw=none},
+ legend cell align = left,
+ xlabel= {iterations},
+ ylabel= {subsequent errors},
+% xmin= 0,
+% xmax= 53,
+% ymin= 0,
+% ymax= 0.0003,
+ grid= both, %major or minor
+ axis line style={-Latex[round]},
+ legend style={
+% anchor=north east,
+% at={(1,1)},
+ font=\tiny
+ },
+ legend entries={$\bigl\|p_l^i - p_l^{i-1}\bigr\|_{L^2(\dom_1)}$,
+ $\bigl\|p_l^i - p_l^{i-1}\bigr\|_{L^2(\dom_2)}$},
+ ] % end of axis options %col sep=comma,
+ \addplot table [x=iteration, y=wetting] {../output/mesh_res30_dt0.004/subdomain1_dt0.004_hmin0.019_lambda70_Lw0.25subsequent_iteration_errors_at_time0.10.csv};%
+ \addplot table [x=iteration, y=wetting] {../output/mesh_res30_dt0.004/subdomain2_dt0.004_hmin0.019_lambda70_Lw0.25subsequent_iteration_errors_at_time0.10.csv};%
+% \addplot table [col sep=comma] {../output/subdomain1_subsequent_iteration_error_for_phase_wetting_at_time0.csv};%
+% % \addplot table [col sep=comma] {../output/subdomain2_subsequent_iteration_error_for_phase_wetting_at_time0.csv};%
+ \end{semilogyaxis}
+ \end{tikzpicture}
+
+ % Generation of Subsequent_errors_t0.00
+ \tikzsetnextfilename{Subsequent_errors_t0.14}
+ \tikzset{external/force remake}
+ \begin{tikzpicture}
+% % \path[draw,dashed,thick] (3.13,0) -- (3.13,5.2);
+ \begin{semilogyaxis}[%
+ width=\textwidth,
+ title={ Subsequent errors for $t = 0.14$, $ h \approx 0.02$, $\tau = 4\cdot 10^{-4}$ },
+ % axis lines=left,
+ % legend style = {draw=none},
+ legend cell align = left,
+ xlabel= {iterations},
+ ylabel= {subsequent errors},
+% xmin= 0,
+% xmax= 53,
+% ymin= 0,
+% ymax= 0.0003,
+ grid= both, %major or minor
+ axis line style={-Latex[round]},
+ legend style={
+% anchor=north east,
+% at={(1,1)},
+ font=\tiny
+ },
+ legend entries={$\bigl\|p_l^i - p_l^{i-1}\bigr\|_{L^2(\dom_1)}$,
+ $\bigl\|p_l^i - p_l^{i-1}\bigr\|_{L^2(\dom_2)}$},
+ ] % end of axis options %col sep=comma,
+ \addplot table [x=iteration, y=wetting] {../output/mesh_res30_dt0.004/subdomain1_dt0.004_hmin0.019_lambda70_Lw0.25subsequent_iteration_errors_at_time0.14.csv};%
+ \addplot table [x=iteration, y=wetting] {../output/mesh_res30_dt0.004/subdomain2_dt0.004_hmin0.019_lambda70_Lw0.25subsequent_iteration_errors_at_time0.14.csv};%
% \addplot table [col sep=comma] {../output/subdomain1_subsequent_iteration_error_for_phase_wetting_at_time0.csv};%
% % \addplot table [col sep=comma] {../output/subdomain2_subsequent_iteration_error_for_phase_wetting_at_time0.csv};%
\end{semilogyaxis}
diff --git a/TP-TP-patch-test-case/TP-TP-2-patch-test.py b/TP-TP-patch-test-case/TP-TP-2-patch-test.py
index 1e9a7d5..6f120ff 100755
--- a/TP-TP-patch-test-case/TP-TP-2-patch-test.py
+++ b/TP-TP-patch-test-case/TP-TP-2-patch-test.py
@@ -87,10 +87,37 @@ isRichards = {
2: False
}
-Tmax = 1
-timestep_size = 4*0.001
-time_interval = [0, Tmax]
+timestep_size = 4*0.001
+number_of_timesteps = 100
+# 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_analise = 5
+starttime = 0
+
+# We want to write out csv files containing subsequent L2 errors of iterations for
+# number_of_timesteps_to_analise many timesteps. We want to do this with
+# analyse_timesteps = np.linspace(starttime, number_of_timesteps, number_of_timesteps_to_analise)
+# see below. Therefor, we need to check if (number_of_timesteps_to_analise-1) is
+# a divisor of number_of_timesteps,
+if number_of_timesteps_to_analise == 0:
+ analyse_timesteps = None
+elif number_of_timesteps_to_analise == 1:
+ analyse_timesteps = [0]
+else:
+ if not number_of_timesteps%(number_of_timesteps_to_analise-1) == 0:
+ # if number_of_timesteps%(number_of_timesteps_to_analise-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_analise = {number_of_timesteps/number_of_timesteps_to_analise}",
+ "is not an integer,\n")
+ new_analysing_interval = int(round(number_of_timesteps/(number_of_timesteps_to_analise-1)))
+ number_of_timesteps = new_analysing_interval*(number_of_timesteps_to_analise-1)
+ print(f"I'm resetting number_of_timesteps = {number_of_timesteps}\n")
+ analyse_timesteps = np.linspace(0, number_of_timesteps, number_of_timesteps_to_analise)
+
+Tmax = number_of_timesteps*timestep_size
+time_interval = [starttime, Tmax]
+print(f"The simulation interval is [{starttime}, {Tmax}]")
viscosity = {#
# subdom_num : viscosity
@@ -116,10 +143,10 @@ L = {#
lambda_param = {#
# subdom_num : lambda parameter for the L-scheme
- 1 : {'wetting' :100,
- 'nonwetting': 100},#
- 2 : {'wetting' :100,
- 'nonwetting': 100}
+ 1 : {'wetting' :140,
+ 'nonwetting': 140},#
+ 2 : {'wetting' :140,
+ 'nonwetting': 140}
}
## relative permeabilty functions on subdomain 1
@@ -333,7 +360,7 @@ simulation.set_parameters(output_dir = "./output/",#
sources = source_expression,#
initial_conditions = initial_condition,#
dirichletBC_Expressions = dirichletBC,#
- exact_solution_case = True,#
+ exact_solution = exact_solution,#
write2file = write_to_file,#
)
@@ -365,7 +392,7 @@ simulation.initialise()
# df.File('./subdomain2_interface_marker.pvd') << simulation.subdomain[2].interface_marker
# df.File('./test_subdomain1_boundary_marker.pvd') << boundary_marker1
# df.File('./test_subdomain2_boundary_marker.pvd') << boundary_marker2
-analyse_timesteps = [0, 100*timestep_size, 1, 2, 3]
+
simulation.run(analyse_timesteps)
# simulation.LDDsolver(time = 0, debug = True, analyse_timestep = True)
# df.info(parameters, True)
--
GitLab