From 508d3d0acda622e057134ad608fb84240d4aa60d Mon Sep 17 00:00:00 2001
From: David Seus <david.seus@ians.uni-stuttgart.de>
Date: Fri, 3 May 2019 13:25:35 +0200
Subject: [PATCH] before attempting to fix calculation of gli terms
---
LDDsimulation/LDDsimulation.py | 53 +++++++----------
LDDsimulation/domainPatch.py | 65 +++++++--------------
RR-2-patch-test-case/RR-2-patch-test.py | 27 +++++----
TP-TP-patch-test-case/TP-TP-2-patch-test.py | 18 +++---
4 files changed, 66 insertions(+), 97 deletions(-)
diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index 1c2bf30..2b432b6 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -239,25 +239,25 @@ class LDDsimulation(object):
raise RuntimeError('Parameters note set!\nYou forgott to run self.set_parameters(**kwds)')
self._add_parameters_to_output_dirname()
self._init_analyse_timesteps()
- df.info(colored("initialise meshes and marker functions ...\n", "orange"))
+ df.info(colored("initialise meshes and marker functions ...\n", "yellow"))
self._init_meshes_and_markers()
- df.info(colored("initialise interfaces ...\n", "orange"))
+ df.info(colored("initialise interfaces ...\n", "yellow"))
self._init_interfaces()
- df.info(colored("initialise subdomains ...\n", "orange"))
+ df.info(colored("initialise subdomains ...\n", "yellow"))
self._init_subdomains()
- df.info(colored("creating xdmf files for each subdomain ...\n", "orange"))
+ df.info(colored("creating xdmf files for each subdomain ...\n", "yellow"))
self._init_solution_files()
# initial values get initialised here to be able to call the solver at different
# timesteps without having to call self.run().
- df.info(colored("create initial values ...", "orange"))
+ df.info(colored("create initial values ...", "yellow"))
self._init_initial_values()
- df.info(colored("create source Expressions ...", "orange"))
+ df.info(colored("create source Expressions ...", "yellow"))
self._init_sources()
# initialise exact solution expression dictionary if we have an exact
# solution case.
self._init_exact_solution_expression()
if self.exact_solution:
- df.info(colored("writing exact solution for all time steps to xdmf files ...", "orange"))
+ df.info(colored("writing exact solution for all time steps to xdmf files ...", "yellow"))
self.write_exact_solution_to_xdmf()
self._initialised = True
@@ -372,7 +372,6 @@ class LDDsimulation(object):
# iterations in, used to analyse the behaviour of the iterates within
# the current time step.
solution_over_iteration_within_timestep = dict()
- subsequent_errors_over_iteration = 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():
@@ -385,20 +384,11 @@ class LDDsimulation(object):
solution_over_iteration_within_timestep.update( #
{ind: SolutionFile(self._mpi_communicator, filename)}
)
- subsequent_errors_over_iteration.update( {ind: dict()} )
- # # we need a dictionries for each phase.
- # for phase in subdomain.has_phases:
- # subsequent_errors_over_iteration[ind].update( {phase: dict()} )
- #
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()
- # update the time of the source terms before starting the iteration.
- # The sources and sinks are the same throughout the iteration.
- # self._eval_sources(time = time,
- # subdomain_index = ind)
self.update_DirichletBC_dictionary(time = time, #
subdomain_index = ind,
)
@@ -409,9 +399,14 @@ class LDDsimulation(object):
# 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})
@@ -437,6 +432,7 @@ class LDDsimulation(object):
self.Lsolver_step(subdomain_index = sd_index,#
debug = debug,
)
+ # bypass error checking for the first iteration.
if iteration > 0:
subsequent_iter_error = dict()
for phase in subdomain.has_phases:
@@ -446,17 +442,17 @@ class LDDsimulation(object):
subsequent_iter_error.update(
{phase: error_calculated}
)
+
# set flag to stop solving for that phase if error_criterion is met.
if subsequent_iter_error[phase] < error_criterion:
self._calc_isdone_for_phase[sd_index].update({phase: True})
if debug:
- print(f"the subsequent error in iteration {iteration} for phase {phase} is {subsequent_iter_error[phase]}")
+ print(f"t = {self.t}: subsequent error in iteration {iteration} for phase {phase} is {subsequent_iter_error[phase]}")
if analyse_timestep:
# 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]}")
+ # with pgfplots
subsequent_error_filename = self.output_dir\
+self.output_filename_parameter_part[sd_index]\
+"subsequent_iteration_errors" +"_at_time"+\
@@ -467,28 +463,21 @@ class LDDsimulation(object):
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.
- print("max(subsequent_iter_error.values()): ", max(subsequent_iter_error.values()))
total_subsequent_iter_error = max(subsequent_iter_error.values())
# check if error criterion has been met
-
+ # print(f" total_subsequent_error in iter {iteration} = {total_subsequent_iter_error }\n")
if debug:
- print(f" total_subsequent_error in iter {iteration} = {total_subsequent_iter_error }\n")
+ # print(f" total_subsequent_error in iter {iteration} = {total_subsequent_iter_error }\n")
print(f"the maximum distance between any to points in a cell mesh.hmax() on subdomain{sd_index} is h={subdomain.mesh.hmax()}")
- print(f"the minimal distance between any to points in a cell mesh.hmin() on subdomain{sd_index} is h={subdomain.mesh.hmin()}")
+ #print(f"the minimal distance between any to points in a cell mesh.hmin() on subdomain{sd_index} is h={subdomain.mesh.hmin()}")
print(f"the error criterion is tol = {error_criterion}")
# wetting_done = self._calc_isdone_for_phase[sd_index]['wetting']
# nonwetting_done = self._calc_isdone_for_phase[sd_index]['nonwetting']
if total_subsequent_iter_error < error_criterion:
- if debug:
- print(f"subdomain{sd_index} reachd error tol = {error_criterion} in iteration {iteration}.")
+ # if debug:
+ print(f"subdomain{sd_index} reached error = {total_subsequent_iter_error} < tol = {error_criterion} after {iteration} iterations.")
self._calc_isdone_for_subdomain[sd_index] = True
# prepare next iteration
diff --git a/LDDsimulation/domainPatch.py b/LDDsimulation/domainPatch.py
index 1ccbbb1..9fff677 100644
--- a/LDDsimulation/domainPatch.py
+++ b/LDDsimulation/domainPatch.py
@@ -21,6 +21,7 @@ import numpy as np
import typing as tp
import ufl as fl
import LDDsimulation
+from termcolor import colored
# Interface = tp.NewType('Interface', tp.any)
# SolutionFile = tp.NewType('LDDsimulation.SolutionFile', tp.object)
#import domainPatch as domainPatch
@@ -184,7 +185,7 @@ class DomainPatch(df.SubDomain):
# Dictionary of FEM function of self.function_space holding the interface
# dofs of the previous iteration of the neighbouring pressure. This is used
# to assemble the gli terms in the method self.calc_gli_term().
- self.neighbouring_interface_pressure = None
+ # self.neighbouring_interface_pressure = None
# valiable holding a solver object of type df.KrylovSolver. It is set the
# first time LDDsimulation.LDDsolver is run.
self.linear_solver = None
@@ -256,8 +257,9 @@ class DomainPatch(df.SubDomain):
porosity = self.porosity
if self.isRichards:
# this assumes that atmospheric pressure has been normalised to 0.
- pc_prev_iter = self.pressure_prev_iter['wetting']
- pc_prev_timestep = self.pressure_prev_timestep['wetting']
+ # and since
+ pc_prev_iter = -self.pressure_prev_iter['wetting']
+ pc_prev_timestep = -self.pressure_prev_timestep['wetting']
if phase == 'nonwetting':
print("CAUTION: You invoked domainPatch.governing_problem() with\n", #
"Parameter phase = 'nonwetting'. But isRichards = True for",
@@ -365,7 +367,7 @@ class DomainPatch(df.SubDomain):
ds = self.ds(interface.marker_value)
# needed for the read_gli_dofs() functions
interface_dofs = self._dof_indices_of_interface[ind]
- # for each interface, this is a list of dofs indices belonging to another
+ # for each interface, the following is a list of dofs indices belonging to another
# interface. The dofs corresponding to these indices must be multiplied
# by 1/2 to to get an average of the gli terms for dofs belonging to
# two different interfaces.
@@ -406,8 +408,8 @@ class DomainPatch(df.SubDomain):
if self.isRichards:
# assuming that we normalised atmospheric pressure to zero,
- # pc = pw in the Richards case.
- pc_prev = p_prev_iter['wetting']
+ # pc = -pw in the Richards case.
+ pc_prev = -p_prev_iter['wetting']
else:
pw_prev = p_prev_iter['wetting']
pnw_prev = p_prev_iter['nonwetting']
@@ -458,8 +460,8 @@ class DomainPatch(df.SubDomain):
# if the neighbour of our subdomain (with index self.subdomain_index)
# assumes a Richards model, we don't have gli terms for the
# nonwetting phase and assemble the terms as zero form.
- # pass
- gli_assembled_tmp[phase] += df.assemble(df.Constant(0)*ds).get_local()
+ pass
+ # gli_assembled_tmp[phase] += df.assemble(df.Constant(0)*ds).get_local()
else:
# in this case the neighbour assumes two-phase flow
# and we need to calculte a gl0 torm for the nonwetting
@@ -498,10 +500,12 @@ class DomainPatch(df.SubDomain):
if debug:
print(f"in iteration {iteration} and neighbour iteration = {neighbour_iter_num}:\n")
print(f"Interface{ind}.gli_term_prev[{subdomain}]['wetting']=\n",interface.gli_term_prev[subdomain]['wetting'])
+
interface.gli_term_prev[subdomain].update(#
{'wetting': interface.gli_term[neighbour]['wetting'].copy()}
)
-
+ # print(colored(f"id(interface.gli_term_prev[{subdomain}][wetting])= {id(interface.gli_term_prev[subdomain]['wetting'])}","red"))
+ # print(colored(f"id(interface.gli_term[{neighbour}][wetting])= {id(interface.gli_term[neighbour]['wetting'])}","red"))
if debug:
print(f"Interface{ind}.gli_term[{neighbour}]['wetting']=\n",interface.gli_term[neighbour]['wetting'])
print(f"Interface{ind}.gli_term_prev[{subdomain}]['wetting']=\n",interface.gli_term_prev[subdomain]['wetting'])
@@ -538,7 +542,8 @@ class DomainPatch(df.SubDomain):
# read neighbouring pressure values from interface and write them to
# the function self.neighbouring_interface_pressure
- interface.write_pressure_dofs(to_function = self.neighbouring_interface_pressure[phase], #
+ pa_prev = df.interpolate(df.Constant(0), self.function_space[phase])
+ interface.write_pressure_dofs(to_function = pa_prev,#self.neighbouring_interface_pressure[phase], #
interface_dofs = interface_dofs[phase],#
dof_to_vert_map = self.dof2vertex[phase],#
local_to_parent_vertex_map = self.parent_mesh_index,#
@@ -548,7 +553,7 @@ class DomainPatch(df.SubDomain):
# use the above to calculate the gli update with
# in the paper p^{i-1}_{3-l}
- pa_prev = self.neighbouring_interface_pressure[phase]
+ # pa_prev = self.neighbouring_interface_pressure[phase]
phi_a = self.testfunction[phase]
gli_pressure_part = df.assemble(-2*dt*Lambda[phase]*pa_prev*phi_a*ds)
gli_p_part = gli_pressure_part.get_local()
@@ -634,44 +639,14 @@ class DomainPatch(df.SubDomain):
self.function_space = dict()
self.trialfunction = dict()
self.testfunction = dict()
- self.neighbouring_interface_pressure = dict()
+ # self.neighbouring_interface_pressure = dict()
for phase in self.has_phases:
self.function_space.update({phase: df.FunctionSpace(self.mesh, 'P', 1)})
self.trialfunction.update({phase: df.TrialFunction(self.function_space[phase])})
self.testfunction.update({phase: df.TestFunction(self.function_space[phase])})
- self.neighbouring_interface_pressure.update(
- {phase: df.interpolate(df.Constant(0), self.function_space[phase])}
- )
- print(f"Python id of self.neighbouring_interface_pressure[{phase}]", id(self.neighbouring_interface_pressure[phase]))
-
- # if self.isRichards:
- # self.function_space = {'wetting': df.FunctionSpace(self.mesh, 'P', 1)}
- # self.trialfunction = {
- # 'wetting': df.TrialFunction(self.function_space['wetting'])
- # }
- # self.testfunction = {
- # 'wetting': df.TestFunction(self.function_space['wetting'])
- # }
- # self.neighbouring_interface_pressure = {
- # 'wetting' : df.interpolate(df.Constant(0), self.function_space['wetting'])
- # }
- # else:
- # self.function_space = {#
- # 'wetting' : df.FunctionSpace(self.mesh, 'P', 1),#
- # 'nonwetting' : df.FunctionSpace(self.mesh, 'P', 1)#
- # }
- # self.trialfunction = {
- # 'wetting' : df.TrialFunction(self.function_space['wetting']),
- # 'nonwetting' : df.TrialFunction(self.function_space['nonwetting'])
- # }
- # self.testfunction = {
- # 'wetting' : df.TestFunction(self.function_space['wetting']),
- # 'nonwetting' : df.TestFunction(self.function_space['nonwetting'])
- # }
- # self.neighbouring_interface_pressure = {#
- # 'wetting' : df.interpolate(df.Constant(0), self.function_space['wetting']),#
- # 'nonwetting' : df.interpolate(df.Constant(0), self.function_space['nonwetting'])
- # }
+ # self.neighbouring_interface_pressure.update(
+ # {phase: df.interpolate(df.Constant(0), self.function_space[phase])}
+ # )
def _init_dof_maps(self) -> None:
""" calculate dof maps and the vertex to parent map.
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 5370b0b..76d6a07 100755
--- a/RR-2-patch-test-case/RR-2-patch-test.py
+++ b/RR-2-patch-test-case/RR-2-patch-test.py
@@ -85,13 +85,13 @@ isRichards = {
}
##################### MESH #####################################################
-mesh_resolution = 10
+mesh_resolution = 40
##################### TIME #####################################################
-timestep_size = 2*0.0001
-number_of_timesteps = 30
+timestep_size = 4*0.0001
+number_of_timesteps = 5000
# 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
+number_of_timesteps_to_analyse = 11
starttime = 0
@@ -115,8 +115,8 @@ L = {#
lambda_param = {#
# subdom_num : lambda parameter for the L-scheme
- 1 : {'wetting' :160},#
- 2 : {'wetting' :160}
+ 1 : {'wetting' :100},#
+ 2 : {'wetting' :100}
}
## relative permeabilty functions on subdomain 1
@@ -147,9 +147,14 @@ relative_permeability = {#
2: subdomain2_rel_perm
}
-def saturation(pressure, subdomain_index):
+# this function needs to be monotonically decreasing in the capillary_pressure.
+# since in the richards case pc = -pw, this becomes as a function of pw a monotonically
+# INCREASING function like in our Richards-Richards paper. However, since we unify
+# the treatment in the code for Richards and two-phase, we need the same requierment
+# for both cases, two-phase and Richards.
+def saturation(capillary_pressure, subdomain_index):
# inverse capillary pressure-saturation-relationship
- return df.conditional(pressure < 0, 1/((1 - pressure)**(1/(subdomain_index + 1))), 1)
+ return df.conditional(capillary_pressure > 0, 1/((1 + capillary_pressure)**(1/(subdomain_index + 1))), 1)
sat_pressure_relationship = {#
1: ft.partial(saturation, subdomain_index = 1),#
@@ -197,7 +202,7 @@ write_to_file = {
}
# initialise LDD simulation class
-simulation = ldd.LDDsimulation(tol = 1E-12)
+simulation = ldd.LDDsimulation(tol = 1E-14)
simulation.set_parameters(output_dir = "./output/",#
subdomain_def_points = subdomain_def_points,#
isRichards = isRichards,#
@@ -213,11 +218,11 @@ simulation.set_parameters(output_dir = "./output/",#
saturation = sat_pressure_relationship,#
starttime = starttime,#
number_of_timesteps = number_of_timesteps,
- number_of_timesteps_to_analise = number_of_timesteps_to_analise,
+ number_of_timesteps_to_analyse = number_of_timesteps_to_analyse,
timestep_size = timestep_size,#
sources = source_expression,#
initial_conditions = initial_condition,#
- dirichletBC_Expressions = dirichletBC,#
+ dirichletBC_expression_strings = dirichletBC,#
exact_solution = exact_solution,#
write2file = write_to_file,#
)
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 7ad8c97..c078425 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 = 100
-timestep_size = 2*0.001
-number_of_timesteps = 1000
+mesh_resolution = 40
+timestep_size = 0.001
+number_of_timesteps = 2000
# 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
@@ -122,9 +122,9 @@ L = {#
lambda_param = {#
# subdom_num : lambda parameter for the L-scheme
1 : {'wetting' :140,
- 'nonwetting': 1000},#
+ 'nonwetting': 2400},#
2 : {'wetting' :140,
- 'nonwetting': 1000}
+ 'nonwetting': 2400}
}
## relative permeabilty functions on subdomain 1
@@ -147,7 +147,7 @@ def rel_perm2w(s):
# relative permeabilty wetting on subdomain2
return s**3
def rel_perm2nw(s):
- # relative permeabilty nonwetting on subdomain2
+ # relative permeabilty nonwetting on subdosym.cos(0.8*t - (0.8*x + 1/7*y))main2
return (1-s)**2
_rel_perm2w = ft.partial(rel_perm2w)
@@ -168,7 +168,7 @@ relative_permeability = {#
# we set alpha = 0, assume m = 1-1/n (see Helmig) and assume that residual saturation is Sw
def saturation(capillary_pressure, n_index, alpha):
# inverse capillary pressure-saturation-relationship
- return df.conditional(capillary_pressure < 0, 1/((1 + (alpha*capillary_pressure)**n_index)**((n_index - 1)/n_index)), 1)
+ return df.conditional(capillary_pressure > 0, 1/((1 + (alpha*capillary_pressure)**n_index)**((n_index - 1)/n_index)), 1)
# S-pc-relation ship. We use the van Genuchten approach, i.e. pc = 1/alpha*(S^{-1/m} -1)^1/n, where
# we set alpha = 0, assume m = 1-1/n (see Helmig) and assume that residual saturation is Sw
@@ -318,7 +318,7 @@ write_to_file = {
# initialise LDD simulation class
-simulation = ldd.LDDsimulation(tol = 1E-12)
+simulation = ldd.LDDsimulation(tol = 1E-14)
simulation.set_parameters(output_dir = "./output/",#
subdomain_def_points = subdomain_def_points,#
isRichards = isRichards,#
@@ -338,7 +338,7 @@ simulation.set_parameters(output_dir = "./output/",#
timestep_size = timestep_size,#
sources = source_expression,#
initial_conditions = initial_condition,#
- dirichletBC_Expressions = dirichletBC,#
+ dirichletBC_expression_strings = dirichletBC,#
exact_solution = exact_solution,#
write2file = write_to_file,#
)
--
GitLab