From d812de874ddeb39aa4d0771add7d4241bd370e05 Mon Sep 17 00:00:00 2001
From: David Seus <david.seus@ians.uni-stuttgart.de>
Date: Tue, 7 May 2019 12:45:32 +0200
Subject: [PATCH] bevor iterationszahl im Solver von 1 losgeht.
---
LDDsimulation/LDDsimulation.py | 166 +++++++++++++-------
LDDsimulation/domainPatch.py | 26 ++-
RR-2-patch-test-case/RR-2-patch-test.py | 8 +-
TP-TP-patch-test-case/TP-TP-2-patch-test.py | 18 +--
4 files changed, 134 insertions(+), 84 deletions(-)
diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index 764ac95..5442a2d 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -162,7 +162,7 @@ class LDDsimulation(object):
'report': False
}
- self.debug_solver = False
+ self.debug_solver = True
# Dictionaries needed by the LDD Solver
# dictionary holding bool variables for each subdomain indicating wether
# minimal error has been achieved, i.e. the calculation can be considered
@@ -358,63 +358,30 @@ class LDDsimulation(object):
""" calulate the solution on all patches for the next timestep
The heart of the simulation. The function LDDsolver implements the LDD
- iteration, that is solves the problem given the initial data at timestep
+ iteration, that solves the problem given the initial data at timestep
time.
+
+ PARAMATERS
+ ---------------------
+ time float, optional time at which to run the solver. This is internally
+ set to self.t if no time is given. This variable
+ is there for development purposes.
+ debug bool, optional debug flag used to toggle the amount of output
+ being displayed.
+ analyse_timestep bool, optional, flag to toggle wether or not to run
+ routines designed to output data within one timestep
+ iteration.
"""
# in case no time is given, use the time of the class.
if time is None:
time = self.t
error_criterion = 1E-8#min(1E-9, subdomain.mesh.hmax()**2*subdomain.timestep_size)
-
- if analyse_timestep:
- # dictionary holding filename strings for each subdomain to store
- # iterations in, used to analyse the behaviour of the iterates within
- # the current time step.
- solution_over_iteration_within_timestep = dict()
- # before the iteration gets started all iteration numbers must be nulled
- # and the self._calc_isdone_for_subdomain dictionary must be set to False
- for ind, subdomain in self.subdomain.items():
- # create xdmf files for writing out iterations if analyse_timestep is
- # given.
- if analyse_timestep:
- 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)}
- )
-
- self._calc_isdone_for_subdomain.update({ind: False})
- # reset all interface[has_interface].current_iteration[ind] = 0, for
- # index has_interface in subdomain.has_interface
- subdomain.null_all_interface_iteration_numbers()
- self.update_DirichletBC_dictionary(time = time, #
- subdomain_index = ind,
- )
- # # this is the beginning of the solver, so iteration must be set
- # # to 0.
- subdomain.iteration_number = 0
- for phase in subdomain.has_phases:
- # evaluate the sources to the current time.
- subdomain.source[phase].t = time
- # set the solution of the previous timestep as new prev_timestep pressure
- # print(colored(f"id(subdomain{subdomain.subdomain_index}.pressure_prev_timestep[{phase}])", "red"), " =\n", f"{id(subdomain.pressure_prev_timestep[phase])}")
- # print(f"subdomain{subdomain.subdomain_index}.pressure_prev_timestep[{phase}].vector()[:]", " =\n", f"{subdomain.pressure_prev_timestep[phase].vector()[:]}")
- # print(f"subdomain{subdomain.subdomain_index}.pressure[{phase}].vector()[:]", " =\n", f"{subdomain.pressure[phase].vector()[:]}")
- subdomain.pressure_prev_timestep[phase].assign(
- subdomain.pressure[phase]
- )
- # print(colored(f"id(subdomain{subdomain.subdomain_index}.pressure[{phase}])", "red"), " =\n", f"{id(subdomain.pressure[phase])}")
- # print(f"subdomain{subdomain.subdomain_index}.pressure_prev_timestep[{phase}].vector()[:]", " =\n", f"{subdomain.pressure_prev_timestep[phase].vector()[:]}")
- # All isdone flags need to be zero before we start iterating
- self._calc_isdone_for_phase[ind].update({phase: False})
-
- subdomain.calc_gl0_term()
-
+ # in case analyse_timestep is None, solution_over_iteration_within_timestep is None
+ solution_over_iteration_within_timestep = self.prepare_LDDsolver(analyse_timestep = analyse_timestep)
### actual iteration starts here
- # gobal stopping criterion for the iteration.
iteration = 0
+ # gobal stopping criterion for the iteration.
all_subdomains_are_done = False
while iteration <= self._max_iter_num and not all_subdomains_are_done:
# we need to loop over the subdomains and solve an L-scheme type step
@@ -426,10 +393,10 @@ class LDDsimulation(object):
if not self._calc_isdone_for_subdomain[sd_index]:
# set subdomain iteration to current iteration
subdomain.iteration_number = iteration
- if debug:
- print(f"On subdomain{sd_index}: Entering iteration {iteration}\n")
- print(f"subdomain{sd_index}.isRichards = {subdomain.isRichards}.")
- print(f"subdomain{sd_index}.has_phases = {subdomain.has_phases}.")
+ # if debug:
+ # print(f"On subdomain{sd_index}: Entering iteration {iteration}\n")
+ # print(f"subdomain{sd_index}.isRichards = {subdomain.isRichards}.")
+ # print(f"subdomain{sd_index}.has_phases = {subdomain.has_phases}.")
# solve the problem on subdomain
self.Lsolver_step(subdomain_index = sd_index,#
debug = debug,
@@ -517,6 +484,97 @@ class LDDsimulation(object):
iteration += 1
# end iteration while loop.
+ def prepare_LDDsolver(self, time: float = None, #
+ debug: bool = False, #
+ analyse_timestep: bool = False) -> analysing_solution_file_dict: tp.Dict[ind, SolutionFile]:
+ """ initialise LDD iteration for current time time.
+
+ This method executes a few prelininary steps before actually starting
+ the L-iteration. Namely,
+ - create solution files for analysing time step if analyse_timestep == True
+ - the dictionries self._calc_isdone_for_subdomain are set to False
+ - Dirichlet boundary condition terms are evaluated to the current time
+ - all subdomain iteration numbers on interfaces are set to 0.
+ - all subdomain iteration numbers subdomain.iteration_number are set
+ to 0.
+ - source terms are evaluated to the current time.
+ - the calculated solution pressure from the previous timestep
+ is assigned to pressure_prev for each subdomain.
+ - The self._calc_isdone_for_phase dictionary for each phase on each
+ subdomain are set to False
+ - gl0 terms are calculated.
+
+ PARAMETERS
+ ------------------------------------
+ time float, optional time at which to run the solver. This is internally
+ set to self.t if no time is given. This variable
+ is there for development purposes.
+ debug bool, optional debug flag used to toggle the amount of output
+ being displayed.
+ analyse_timestep bool, optional, flag to toggle wether or not to run
+ routines designed to output data within one timestep
+ iteration.
+
+ OUTPUT
+ ------------------------------------
+ analysing_solution_file_dict tp.Dict[ind, SolutionFile] dictionary
+ containing the solution file for analysing the
+ timestep if analyse_timestep == True, else None.
+ """
+ # in case no time is given, use the time of the class.
+ if time is None:
+ time = self.t
+
+ if analyse_timestep:
+ # dictionary holding filename strings for each subdomain to store
+ # iterations in, used to analyse the behaviour of the iterates within
+ # the current time step.
+ solution_over_iteration_within_timestep = dict()
+ else:
+ solution_over_iteration_within_timestep = None
+ # before the iteration gets started all iteration numbers must be nulled
+ # and the self._calc_isdone_for_subdomain dictionary must be set to False
+ for ind, subdomain in self.subdomain.items():
+ # create xdmf files for writing out iterations if analyse_timestep is
+ # given.
+ if analyse_timestep:
+ 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)}
+ )
+
+ self._calc_isdone_for_subdomain.update({ind: False})
+ # reset all interface[has_interface].current_iteration[ind] = 0, for
+ # index has_interface in subdomain.has_interface
+ subdomain.null_all_interface_iteration_numbers()
+ self.update_DirichletBC_dictionary(time = time, #
+ subdomain_index = ind,
+ )
+ # # this is the beginning of the solver, so iteration must be set
+ # # to 0.
+ subdomain.iteration_number = 0
+ for phase in subdomain.has_phases:
+ # evaluate the sources to the current time.
+ subdomain.source[phase].t = time
+ # set the solution of the previous timestep as new prev_timestep pressure
+ if debug:
+ print(colored(f"id(subdomain{subdomain.subdomain_index}.pressure_prev_timestep[{phase}])", "red"), " =\n", f"{id(subdomain.pressure_prev_timestep[phase])}")
+ print(f"subdomain{subdomain.subdomain_index}.pressure_prev_timestep[{phase}].vector()[:]", " =\n", f"{subdomain.pressure_prev_timestep[phase].vector()[:]}")
+ print(f"subdomain{subdomain.subdomain_index}.pressure[{phase}].vector()[:]", " =\n", f"{subdomain.pressure[phase].vector()[:]}")
+ subdomain.pressure_prev_timestep[phase].assign(
+ subdomain.pressure[phase]
+ )
+ if debug:
+ print(colored(f"id(subdomain{subdomain.subdomain_index}.pressure[{phase}])", "red"), " =\n", f"{id(subdomain.pressure[phase])}")
+ print(f"subdomain{subdomain.subdomain_index}.pressure_prev_timestep[{phase}].vector()[:]", " =\n", f"{subdomain.pressure_prev_timestep[phase].vector()[:]}")
+ # All isdone flags need to be zero before we start iterating
+ self._calc_isdone_for_phase[ind].update({phase: False})
+
+ subdomain.calc_gl0_term()
+
+ return solution_over_iteration_within_timestep
def Lsolver_step(self, subdomain_index: int,#
debug: bool = False) -> None:
@@ -796,7 +854,7 @@ class LDDsimulation(object):
def _init_interfaces(self, interface_def_points: tp.List[tp.List[df.Point]] = None,#
adjacent_subdomains: tp.List[np.ndarray] = None) -> None:
- """ generate interfaces and interface markerfunctions for suddomains and
+ """ generate interfaces and interface marker functions for suddomains and
set the internal lists used by the LDDsimulation class.
This initialises a list of interfaces and global interface marker functions
diff --git a/LDDsimulation/domainPatch.py b/LDDsimulation/domainPatch.py
index 3c50a01..76832d4 100644
--- a/LDDsimulation/domainPatch.py
+++ b/LDDsimulation/domainPatch.py
@@ -368,28 +368,24 @@ class DomainPatch(df.SubDomain):
f"dof2global_vertex_map:\n{self.parent_mesh_index[self.dof2vertex['wetting'][interface_dofs['wetting']]]}\n"
f"and dofs common with other interfaces:\n{dofs_in_common_with_other_interfaces['wetting']}")
- n = df.FacetNormal(self.mesh)
- pa_prev_timestep = self.pressure_prev_timestep
- k = self.relative_permeability
- S = self.saturation
- mu = self.viscosity
-
if self.isRichards:
# assuming that we normalised atmospheric pressure to zero,
# pc = -pw in the Richards case.
- pc_prev = -pa_prev_timestep['wetting']
+ pc_prev = -self.pressure_prev_timestep['wetting']
else:
- pw_prev = pa_prev_timestep['wetting']
- pnw_prev = pa_prev_timestep['nonwetting']
+ pw_prev = self.pressure_prev_timestep['wetting']
+ pnw_prev = self.pressure_prev_timestep['nonwetting']
pc_prev = pnw_prev - pw_prev
+ n = df.FacetNormal(self.mesh)
+ S = self.saturation
for phase in self.has_phases:
# viscosity
- mu_a = mu[phase]
+ mu_a = self.viscosity[phase]
# relative permeability
- ka = k[phase]
+ ka = self.relative_permeability[phase]
# previous timestep pressure
- pa_prev = pa_prev_timestep[phase]
+ pa_prev = self.pressure_prev_timestep[phase]
# testfunction
phi_a = self.testfunction[phase]
if phase == 'nonwetting' and interface.neighbour_isRichards[subdomain]:
@@ -443,9 +439,6 @@ class DomainPatch(df.SubDomain):
# save the result of the calculation to the subdomain and to the interface
# dictionaries.
for phase in self.has_phases:
- # self.gli_assembled.update(
- # {phase: gli_assembled_tmp[phase].copy()}
- # )
for interf_ind in self.has_interface:
# self._calc_gli_term_on(interface, iteration)
interface = self.interface[interf_ind]
@@ -516,8 +509,7 @@ class DomainPatch(df.SubDomain):
# in case phase == 'nonwetting' only proceed if the neighbouring
# subdomain does not assume a Richards model.
if (phase == 'wetting') or (not interface.neighbour_isRichards[subdomain]):
- # read gli_prev term from the neighbour and write it to
- # the gli_assembled array of the subdomain
+ # read gli_prev term from the neighbour.
# to be able to sum the gli terms of all interfaces together,
# we need a dummy vector to hold the dofs read from the interface.
gli_prev_of_neighbour = np.zeros(gli_assembled_tmp[phase].size, dtype = float)
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 728b86c..4cc538c 100755
--- a/RR-2-patch-test-case/RR-2-patch-test.py
+++ b/RR-2-patch-test-case/RR-2-patch-test.py
@@ -85,10 +85,10 @@ isRichards = {
}
##################### MESH #####################################################
-mesh_resolution = 6
+mesh_resolution = 100
##################### TIME #####################################################
timestep_size = 4*0.0001
-number_of_timesteps = 200
+number_of_timesteps = 3000
# 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 = 11
@@ -115,8 +115,8 @@ L = {#
lambda_param = {#
# subdom_num : lambda parameter for the L-scheme
- 1 : {'wetting' :4},#
- 2 : {'wetting' :4}
+ 1 : {'wetting' :100},#
+ 2 : {'wetting' :100}
}
## relative permeabilty functions on subdomain 1
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 c078425..daf7b63 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
@@ -89,9 +89,9 @@ isRichards = {
############ GRID ########################ΓΌ
-mesh_resolution = 40
-timestep_size = 0.001
-number_of_timesteps = 2000
+mesh_resolution = 20
+timestep_size = 0.0001
+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_analyse = 11
@@ -113,18 +113,18 @@ porosity = {#
L = {#
# subdom_num : subdomain L for L-scheme
- 1 : {'wetting' :0.3,
+ 1 : {'wetting' :0.4,
'nonwetting': 0.8},#
- 2 : {'wetting' :0.3,
+ 2 : {'wetting' :0.4,
'nonwetting': 0.8}
}
lambda_param = {#
# subdom_num : lambda parameter for the L-scheme
- 1 : {'wetting' :140,
- 'nonwetting': 2400},#
- 2 : {'wetting' :140,
- 'nonwetting': 2400}
+ 1 : {'wetting' :10,
+ 'nonwetting': 100},#
+ 2 : {'wetting' :10,
+ 'nonwetting': 100}
}
## relative permeabilty functions on subdomain 1
--
GitLab