From b4aa88dfc69cdd0e1bb7c0d9e77cdf65e128a381 Mon Sep 17 00:00:00 2001
From: David Seus <david.seus@ians.uni-stuttgart.de>
Date: Fri, 16 Aug 2019 14:32:08 +0200
Subject: [PATCH] fix wrong assignment of gli_dofs on interface facets
---
.gitignore | 2 +
LDDsimulation/LDDsimulation.py | 527 +++++++++++++++++-------
LDDsimulation/boundary_and_interface.py | 39 +-
LDDsimulation/domainPatch.py | 388 +++++++++++------
LDDsimulation/helpers.py | 283 +++++--------
LDDsimulation/solutionFile.py | 2 +-
6 files changed, 781 insertions(+), 460 deletions(-)
diff --git a/.gitignore b/.gitignore
index d66cabb..eb49951 100644
--- a/.gitignore
+++ b/.gitignore
@@ -57,6 +57,8 @@
*.h5
*.dep
*.csv
+*.ipynb
+core
# Ignoriere Bilder und Graphiken sowie Videos und Musik
*.png
diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index 71b503d..455b2e7 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -1,10 +1,12 @@
-"""
-LDD simulation class
+"""LDD simulation class.
+
+Main class of the LDD simulation. The heart and soul of the code.
"""
import os
import dolfin as df
import numpy as np
+import numpy.linalg as la
import typing as tp
import mshr
import domainPatch as dp
@@ -19,9 +21,8 @@ from termcolor import colored
from solutionFile import SolutionFile
-
class LDDsimulation(object):
- """ Main Class of LDD simulation
+ """Main Class of LDD simulation.
LDDsimulation is the main class for the LDD simulation. It contains methods for
setting up and running the LDD simulation for different meshes and domains.
@@ -54,14 +55,16 @@ class LDDsimulation(object):
# # To be able to run DG in parallel
# df.parameters["ghost_mode"] = "shared_facet"
# df.parameters["ghost_mode"] = "none"
+ df.parameters["mesh_partitioner"] = "ParMETIS"
# Form compiler options
df.parameters["form_compiler"]["optimize"] = True
df.parameters["form_compiler"]["cpp_optimize"] = True
# df.parameters['std_out_all_processes'] = False
+ # df.set_log_level(df.LogLevel.WARNING)
df.set_log_level(df.LogLevel.INFO)
- # df.set_log_level(df.LogLevel.INFO)
+ # df.set_log_level(df.LogLevel.DEBUG)
# CRITICAL = 50, // errors that may lead to data corruption and suchlike
# ERROR = 40, // things that go boom
@@ -69,7 +72,7 @@ class LDDsimulation(object):
# INFO = 20, // information of general interest
# PROGRESS = 16, // what's happening (broadly)
# TRACE = 13, // what's happening (in detail)
- # DBG = 10 // sundry
+ # DEBUG = 10, // sundry
hlp.print_once("Set default parameter...\n")
self.restart_file = None
@@ -92,7 +95,7 @@ class LDDsimulation(object):
# global marker for subdomains and interfaces. Is set by self.init_meshes_and_markers()
self.domain_marker = None
# List of interfaces of class Interface. This is populated by self.init_interfaces
- self.interface = []
+ self.interface = None
# The global interface marker. Gets initialised by self.init_interfaces().
self.interface_marker = None
# List of list of indices indicating subdomains adjacent to interfaces.
@@ -103,7 +106,7 @@ class LDDsimulation(object):
## Private variables
# maximal number of L-iterations that the LDD solver uses.
- self._max_iter_num = 100
+ self._max_iter_num = 500
# TODO rewrite this with regard to the mesh sizes
# self.calc_tol = self.tol
# list of timesteps that get analyed. Gets initiated by self._init_analyse_timesteps
@@ -145,8 +148,8 @@ class LDDsimulation(object):
# subdomain. This gets populated by self._init_solution_files()
self.solution_file = None
########## SOLVER DEFINITION AND PARAMETERS ########
- # print("\nLinear Algebra Backends:")
- # df.list_linear_algebra_backends()
+ print("\nLinear Algebra Backends:")
+ df.list_linear_algebra_backends()
# # print("\nLinear Solver Methods:")
df.list_linear_solver_methods()
print("\n")
@@ -155,24 +158,37 @@ class LDDsimulation(object):
# # print("\nPeconditioners for Krylov Solvers:")
df.list_krylov_solver_preconditioners()
- # df.LinearSolver_default_parameters()
+
df.info(df.parameters, True)
- self.solver_type_is_Kryov = False
+ self.solver_type_is_Krylov = True
### Define the linear solver to be used.
- self.solver = 'superlu' #'gmres'#'bicgstab' # biconjugate gradient stabilized method
- self.preconditioner = 'default'#jacobi#'hypre_amg' #'ilu'#'hypre_amg' # incomplete LU factorization
- # dictionary of solver parametrs. This is passed to self._init_subdomains,
- # where for each subdomain a sovler object of type self.solver is created
- # with these parameters.
- self.solver_parameters = {
- 'nonzero_initial_guess': True,
- 'absolute_tolerance': 1E-14,
- 'relative_tolerance': 1E-12,
- 'maximum_iterations': 1000,
- 'report': False
- }
+ if self.solver_type_is_Krylov:
+ self.solver = 'bicgstab'#'superlu' #'gmres'#'bicgstab' # biconjugate gradient stabilized method
+ self.preconditioner = 'jacobi' #'hypre_amg' 'default' #'ilu'#'hypre_amg' # incomplete LU factorization
+ # dictionary of solver parametrs. This is passed to self._init_subdomains,
+ # where for each subdomain a sovler object of type self.solver is created
+ # with these parameters.
+ self.solver_parameters = {
+ 'nonzero_initial_guess': True,
+ 'absolute_tolerance': 1E-14,
+ 'relative_tolerance': 1E-12,
+ 'maximum_iterations': 1000,
+ 'report': False,
+ 'monitor_convergence': False
+ }
+ else:
+ self.solver = 'superlu' #'gmres'#'bicgstab' # biconjugate gradient stabilized method
+ self.preconditioner = 'default' #'hypre_amg' 'default' #'ilu'#'hypre_amg' # incomplete LU factorization
+ # dictionary of solver parametrs. This is passed to self._init_subdomains,
+ # where for each subdomain a sovler object of type self.solver is created
+ # with these parameters.
+ self.solver_parameters = {
+ 'report' : True,
+ 'symmetric': False,
+ 'verbose': True
+ }
self.debug_solver = debug
self.LDDsolver_tol = LDDsolver_tol
@@ -182,6 +198,7 @@ class LDDsimulation(object):
# done or not.
def set_parameters(self,#
+ use_case: str,
output_dir: str,#
subdomain_def_points: tp.List[tp.List[df.Point]],#
# this is actually an array of bools
@@ -212,6 +229,9 @@ class LDDsimulation(object):
)-> None:
""" set parameters of an instance of class LDDsimulation"""
+ self.use_case = use_case
+ if include_gravity:
+ self.use_case = self.use_case + "-with-gravity"
self.output_dir = output_dir
self.subdomain_def_points = subdomain_def_points
self.isRichards = isRichards
@@ -277,7 +297,10 @@ class LDDsimulation(object):
self.write_exact_solution_to_xdmf()
self._initialised = True
- def run(self, analyse_solver_at_timesteps: tp.List[float] = None):
+
+ def run(self,
+ analyse_solver_at_timesteps: tp.List[float] = None,
+ analyse_condition: bool = False):
""" run the simulation
This method implements the time stepping, calls the acutal solvers and writes
@@ -309,24 +332,31 @@ class LDDsimulation(object):
df.info(colored("Start time stepping","green"))
# write initial values to file.
- self.timestep_num = int(0)
+ 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:
+ print("LDD simulation use case: {}".format(self.use_case))
print(f"entering timestep ",
"{}".format(self.timestep_num),
"at time t = "+ "{number:.{digits}f}".format(number = self.t, digits = 4))
# check if the solver should beanalised or not.
if analyse_solver_at_timesteps is not None and np.isin(self.timestep_num, analyse_solver_at_timesteps, assume_unique=True):
print("analyising timestep ...")
- self.LDDsolver(debug = self.debug_solver, analyse_timestep = True)
+ self.LDDsolver(debug=self.debug_solver,
+ analyse_timestep=True,
+ analyse_condition=analyse_condition)
else:
- self.LDDsolver(debug = self.debug_solver, analyse_timestep = False)
+ self.LDDsolver(debug=self.debug_solver,
+ analyse_timestep=False,
+ analyse_condition=False)
self.t += self.timestep_size
- self.timestep_num += 1
# write solutions to files.
for subdom_ind, subdomain in self.subdomain.items():
+ # for phase in subdomain.has_phases:
+ # print(f"calculated pressure of {phase}phase:")
+ # print(subdomain.pressure[phase].vector()[:])
subdomain.write_solution_to_xdmf(#
file = self.solution_file[subdom_ind], #
time = self.t,#
@@ -336,11 +366,29 @@ class LDDsimulation(object):
# absolute error.
if subdomain.pressure_exact is not None:
relative_L2_errornorm = dict()
+ pressure_exact = dict()
for phase in subdomain.has_phases:
pa_exact = subdomain.pressure_exact[phase]
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)
+ pressure_exact.update(
+ {phase: df.interpolate(pa_exact, subdomain.function_space["pressure"][phase])}
+ )
+ absolute_difference = df.Function(subdomain.function_space["pressure"][phase])
+ pressure_difference = pressure_exact[phase].vector()[:] - subdomain.pressure[phase].vector()[:]
+ abs_diff_tmp = np.fabs(pressure_difference)
+ absolute_difference.vector()[:] = abs_diff_tmp
+ if phase == "wetting":
+ data_set_name="abs_diff_w"
+ else:
+ data_set_name="abs_diff_nw"
+ self.write_function_to_xdmf(#
+ function=absolute_difference,
+ file=self.solution_file[subdom_ind], #
+ time=self.t,#
+ data_set_name=data_set_name
+ )
if norm_pa_exact > self.tol:
relative_L2_errornorm.update(
{phase: error_calculated/norm_pa_exact}
@@ -357,9 +405,10 @@ class LDDsimulation(object):
errors = relative_L2_errornorm,
)
+ self.timestep_num += 1
-
- df.info("Finished")
+ print("LDD simulation use case: {}".format(self.use_case))
+ df.info("Finished calculation")
df.info("Elapsed time:"+str(timer.elapsed()[0]))
timer.stop()
# r_cons.close()
@@ -370,12 +419,15 @@ class LDDsimulation(object):
# with open(self.output_dir+"timings.txt", "w") as timefile:
# timefile.write(df.timings)
- df.dump_timings_to_xml(self.output_dir+"timings.xml",df.TimingClear.keep)
-
+ df.dump_timings_to_xml(self.output_dir+"timings.xml", df.TimingClear.keep)
+ # df.info(colored("start post processing calculations ...\n", "yellow"))
+ # self.post_processing()
+ # df.info(colored("finished post processing calculations \nAll right. I'm Done.", "green"))
def LDDsolver(self, time: float = None, #
debug: bool = False, #
- analyse_timestep: bool = False):
+ analyse_timestep: bool = False,
+ analyse_condition: bool = False):
""" calulate the solution on all patches for the next timestep
The heart of the simulation. The function LDDsolver implements the LDD
@@ -408,7 +460,7 @@ class LDDsimulation(object):
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:
+ while iteration < self._max_iter_num and not all_subdomains_are_done:
iteration += 1
# we need to loop over the subdomains and solve an L-scheme type step
# on each subdomain. here it should be possible to parallelise in
@@ -417,44 +469,62 @@ class LDDsimulation(object):
# set subdomain iteration to current iteration
subdomain.iteration_number = iteration
# solve the problem on subdomain
- self.Lsolver_step(subdomain_index = sd_index,#
- debug = debug,
+ condition_nums = self.Lsolver_step(
+ subdomain_index=sd_index,#
+ debug=debug,
+ analyse_condition=analyse_condition
)
# bypass subsequent error calculation for the first iteration, becauses
# it is zero anyways.
- if iteration > 1:
- subsequent_iter_error = dict()
- for phase in subdomain.has_phases:
- error = df.Function(subdomain.function_space["pressure"][phase])
- error.assign(subdomain.pressure[phase] - subdomain.pressure_prev_iter[phase])
- error_calculated = df.norm(error, 'L2')
- subsequent_iter_error.update(
- {phase: error_calculated}
- )
- if debug:
- print(f"time = {time}: subsequent error on subdomain {sd_index} for {phase} phase in iteration {iteration} = ", subsequent_iter_error[phase])
+ #if iteration >= 1:
+ subsequent_iter_error = dict()
+ if debug:
+ print("LDD simulation use case: {}".format(self.use_case))
+ for phase in subdomain.has_phases:
+ error = df.Function(subdomain.function_space["pressure"][phase])
+ error.assign(subdomain.pressure[phase] - subdomain.pressure_prev_iter[phase])
+ error_calculated = df.norm(error, 'L2')
+ subsequent_iter_error.update(
+ {phase: error_calculated}
+ )
+ if debug:
+ print(f"time = {time} (at timestep {self.timestep_num}): subsequent error on subdomain {sd_index} for {phase} phase in iteration {iteration} = ", subsequent_iter_error[phase])
- # calculate the maximum over phase of subsequent errors for the
- # stopping criterion.
- max_subsequent_error[sd_index-1] = max(subsequent_iter_error.values())
+ # calculate the maximum over phase of subsequent errors for the
+ # stopping criterion.
+ max_subsequent_error[sd_index-1] = max(subsequent_iter_error.values())
- if analyse_timestep:
- # save the calculated L2 errors to csv file for plotting
- # with pgfplots
- subsequent_error_filename = self.output_dir\
+ if analyse_timestep:
+ # save the calculated L2 errors to csv file for plotting
+ # with pgfplots
+ subsequent_error_filename = self.output_dir\
+ +self.output_filename_parameter_part[sd_index]\
+ +"subsequent_iteration_errors" +"_at_time"+\
+ "{number}".format(number=self.timestep_num) +".csv" #"{number:.{digits}f}".format(number=time, digits=4)
+ self.write_subsequent_errors_to_csv(
+ filename = subsequent_error_filename, #
+ subdomain_index = sd_index,
+ errors = subsequent_iter_error
+ )
+ if analyse_condition:
+ # save the calculated condition numbers of the assembled
+ # matrices a separate file for monitoring
+ condition_number_filename = self.output_dir\
+self.output_filename_parameter_part[sd_index]\
- +"subsequent_iteration_errors" +"_at_time"+\
- "{number:.{digits}f}".format(number=time, digits=4) +".csv"
- self.write_subsequent_errors_to_csv(
- filename = subsequent_error_filename, #
+ +"condition_numbers" +"_at_timestep"+\
+ "{number}".format(number=self.timestep_num) +".csv"
+ self.write_condition_numbers_to_csv(
+ filename = condition_number_filename, #
subdomain_index = sd_index,
- errors = subsequent_iter_error,
- time = time)
+ condition_number = condition_nums
+ )
# prepare next iteration
# write the newly calculated solution to the inteface dictionaries
# of the subdomain for communication.
- subdomain.write_pressure_to_interfaces()
+ if self.interface_def_points is not None:
+ subdomain.write_pressure_to_interfaces()
+
if analyse_timestep:
subdomain.write_solution_to_xdmf(#
file = solution_over_iteration_within_timestep[sd_index], #
@@ -469,22 +539,23 @@ class LDDsimulation(object):
# end loop over subdomains
##### stopping criterion for the solver.
# only check if error criterion has been met after at least one iteration.
- if iteration > 1:
- if debug:
- # 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 error criterion is tol = {error_criterion}")
-
- total_subsequent_error = np.amax(max_subsequent_error)
- if debug:
- print(f"time = {time}: max subsequent error of both subdomain in iteration {iteration} = ", total_subsequent_error)
- if total_subsequent_error < error_criterion:
- # if debug:
- print(f"all phases on all subdomains reached a subsequent error = {total_subsequent_error} < tol = {error_criterion} after {iteration} iterations.")
- all_subdomains_are_done = True
+ #if iteration >= 1:
+ if debug:
+ # 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 error criterion is tol = {error_criterion}")
+
+ total_subsequent_error = np.amax(max_subsequent_error)
+ if debug:
+ print(f"time = {time} (at timestep {self.timestep_num}):: max subsequent error of both subdomain in iteration {iteration} = ", total_subsequent_error)
+ if total_subsequent_error < error_criterion:
+ # if debug:
+ print(f"all phases on all subdomains reached a subsequent error = {total_subsequent_error} < tol = {error_criterion} after {iteration} iterations.")
+ all_subdomains_are_done = True
# end iteration while loop.
+
def prepare_LDDsolver(self, time: float = None, #
debug: bool = False, #
analyse_timestep: bool = False) -> tp.Dict[int, tp.Type[SolutionFile]]:
@@ -545,7 +616,9 @@ class LDDsimulation(object):
# reset all interface[has_interface].current_iteration[ind] = 0, for
# index has_interface in subdomain.has_interface
- subdomain.null_all_interface_iteration_numbers()
+ if self.interface_def_points is not None:
+ subdomain.null_all_interface_iteration_numbers()
+
if subdomain.outer_boundary is not None:
self.update_DirichletBC_dictionary(time = time, #
subdomain_index = ind,
@@ -572,14 +645,17 @@ class LDDsimulation(object):
# method works properly. The current pressure dictionaries should already
# be populated with the current solution, sinces the solver writes the pressure
# to interfaces after each iteration.
- subdomain.write_pressure_to_interfaces()
- subdomain.write_pressure_to_interfaces(previous_iter=True)
- subdomain.calc_gl0_term()
+ if self.interface_def_points is not None:
+ subdomain.write_pressure_to_interfaces()
+ subdomain.write_pressure_to_interfaces(previous_iter=True)
+ subdomain.calc_gl0_term()
return solution_over_iteration_within_timestep
+
def Lsolver_step(self, subdomain_index: int,#
- debug: bool = False) -> None:
+ analyse_condition: bool = False,
+ debug: bool = False) -> tp.Dict[str, float]:
""" L-scheme solver iteration step for an object of class subdomain
Lsolver_step implements L-scheme solver iteration step for an subdomain
@@ -596,7 +672,9 @@ class LDDsimulation(object):
"""
subdomain = self.subdomain[subdomain_index]
iteration = subdomain.iteration_number
-
+ condition_number = None
+ if analyse_condition:
+ condition_number = dict()
for phase in subdomain.has_phases:
# extract L-scheme form and rhs (without gli term) from subdomain.
governing_problem = subdomain.governing_problem(phase = phase)
@@ -606,11 +684,10 @@ class LDDsimulation(object):
form_assembled = df.assemble(form)
rhs_assembled = df.assemble(rhs)
- if debug and subdomain.mesh.num_cells() < 36:
- # print(f"phase = {phase}: form_assembled:\n", form_assembled.array())
- # print(f"phase = {phase}: rhs:\n", rhs_assembled.get_local())
- print(f"phase = {phase}: gli_term:\n", gli_term_assembled)
- print(f"phase = {phase}: rhs_assembled:\n", rhs_assembled.get_local())
+ if analyse_condition:
+ condition_number.update(
+ {phase: la.cond(form_assembled.array())}
+ )
# apply outer Dirichlet boundary conditions if present.
if subdomain.outer_boundary is not None:
@@ -633,6 +710,8 @@ class LDDsimulation(object):
if debug and subdomain.mesh.num_cells() < 36:
print("\npressure after solver:\n", subdomain.pressure[phase].vector().get_local())
+ return condition_number
+
def _init_solution_files(self):
""" set up solution files for saving the solution of the LDD scheme for
@@ -641,6 +720,7 @@ class LDDsimulation(object):
# set up solution file names for each subdomain.
self.output_filename_parameter_part = dict()
self.solution_file = dict()
+ self.postprocessed_solution_file = dict()
for subdom_ind, subdomain in self.subdomain.items():
tau = subdomain.timestep_size
L = subdomain.L
@@ -653,7 +733,6 @@ class LDDsimulation(object):
+ "_gravitiy_" + "{}".format(self.include_gravity)
name3 = "_lambda" + "{}".format(Lam) + "_Lw" + "{}".format(L["wetting"])
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']
@@ -665,7 +744,9 @@ class LDDsimulation(object):
+"_lambda_nw" + "{}".format(Lam_nw)\
+"_Lw" + "{}".format(L["wetting"])\
+"_Lnw" + "{}".format(L["nonwetting"])
- filename = self.output_dir + filename_param_part + "_solution" +".xdmf"
+
+ filename = self.output_dir + filename_param_part + "_solution" +".xdmf"
+ post_filename = self.output_dir + filename_param_part + "_solution_post" +".xdmf"
# create solution files and populate dictionary.
self.output_filename_parameter_part.update(
{subdom_ind: filename_param_part}
@@ -673,6 +754,10 @@ class LDDsimulation(object):
self.solution_file.update( #
{subdom_ind: SolutionFile(self._mpi_communicator, filename)}
)
+ # self.postprocessed_solution_file.update(
+ # {subdom_ind: SolutionFile(self._mpi_communicator, post_filename)}
+ # )
+
def write_exact_solution_to_xdmf(self):
"""
@@ -682,14 +767,26 @@ class LDDsimulation(object):
for subdom_ind, subdomain in self.subdomain.items():
t = copy.copy(self.starttime)
file = self.solution_file[subdom_ind]
- for timestep in range(0,self.number_of_timesteps+1):
+
+ for timestep in range(self.number_of_timesteps+1):
+ exact_pressure = dict()
for phase in subdomain.has_phases:
pa_exact = subdomain.pressure_exact[phase]
pa_exact.t = t
- tmp_exact_pressure = df.interpolate(pa_exact, subdomain.function_space["pressure"][phase])
- tmp_exact_pressure.rename("exact_pressure_"+"{}".format(phase), "exactpressure_"+"{}".format(phase))
- file.write(tmp_exact_pressure, t)
-
+ exact_pressure.update(
+ {phase: df.interpolate(pa_exact, subdomain.function_space["pressure"][phase])}
+ )
+ exact_pressure[phase].rename("exact_pressure_"+"{}".format(phase), "exact_pressure_"+"{}".format(phase))
+ file.write(exact_pressure[phase], t)
+
+ exact_capillary_pressure = df.Function(subdomain.function_space["pressure"]['wetting'])
+ if self.isRichards:
+ exact_capillary_pressure.assign(-exact_pressure['wetting'])
+ else:
+ pc_temp = exact_pressure["nonwetting"].vector()[:] - exact_pressure["wetting"].vector()[:]
+ exact_capillary_pressure.vector()[:] = pc_temp
+ exact_capillary_pressure.rename("pc_exact", "pc_exact")
+ file.write(exact_capillary_pressure, t)
t += self.timestep_size
@@ -697,12 +794,8 @@ class LDDsimulation(object):
filename: str, #
subdomain_index: int,#
errors: tp.Dict[str, float],#
- time: float = None) -> None:
+ ) -> 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]
- if time is None:
- time = self.t
subdomain = self.subdomain[subdomain_index]
iteration = subdomain.iteration_number
@@ -727,6 +820,37 @@ class LDDsimulation(object):
else:
self.errors.to_csv(filename, sep='\t', encoding='utf-8', index=False)
+
+ def write_condition_numbers_to_csv(self,#
+ filename: str, #
+ subdomain_index: int,#
+ condition_number: tp.Dict[str, float],#
+ ) -> None:
+ """ write subsequent errors to csv file for plotting with pgfplots"""
+ subdomain = self.subdomain[subdomain_index]
+ iteration = subdomain.iteration_number
+
+ if subdomain.isRichards:
+ condition = pd.DataFrame({ "iteration": iteration,
+ #"th. initial mass": df.assemble(rho*df.dx), \
+ #"energy": self.total_energy(), \
+ "wetting": condition_number["wetting"]}, index=[iteration])
+
+ else:
+ condition = pd.DataFrame({"iteration":iteration, #\
+ #"th. initial mass": df.assemble(rho*df.dx), \
+ #"energy": self.total_energy(), \
+ "wetting": condition_number["wetting"],
+ "nonwetting": condition_number["nonwetting"]}, index=[iteration])
+ # if self._rank == 0:
+ # check if file exists
+ if os.path.isfile(filename) == True:
+ with open(filename, 'a') as f:
+ condition.to_csv(f, header=False, sep='\t', encoding='utf-8', index=False)
+ else:
+ condition.to_csv(filename, sep='\t', encoding='utf-8', index=False)
+
+
def write_errornorms_to_csv(self,#
filename: str, #
subdomain_index: int,#
@@ -761,6 +885,28 @@ class LDDsimulation(object):
else:
self.errors.to_csv(filename, sep='\t', encoding='utf-8', index=False)
+
+ def write_function_to_xdmf(self,
+ function: df.Function,
+ file: tp.Type[SolutionFile], #
+ data_set_name: str,
+ time: float = None,#
+ ):
+ """
+ Write function to disk file file. It is expected, that file is a df.XDMFFile
+
+ PARAMETERS
+ function df.Function dolfin function containing the
+ information to write to the file.
+ file str File name of df.XDMFFile to write to
+ data_set_name str name of the data set
+ time float time at which to write the solution.
+
+ """
+ function.rename("{}".format(data_set_name), "{}".format(data_set_name))
+ file.write(function, time)
+
+
## Private methods
def _init_meshes_and_markers(self, subdomain_def_points: tp.List[tp.List[df.Point]] = None,#
mesh_resolution: float = None, #
@@ -806,29 +952,39 @@ class LDDsimulation(object):
### generate global mesh and subdomain meshes. Also count number of subdomains.
# define the global simulation domain polygon first.
domain = mshr.Polygon(subdomain_def_points[0])
- for num, subdomain_vertices in enumerate(subdomain_def_points[1:], 1):
- subdomain = mshr.Polygon(subdomain_vertices)
- domain.set_subdomain(num, subdomain)
- # count number of subdomains for further for loops
- self._number_of_subdomains = num
+ if self.interface_def_points is None:
+ # in this case we don't have any interface and no domain decomposi-
+ # tion going on.
+ self._number_of_subdomains = 1
+ else:
+ for num, subdomain_vertices in enumerate(subdomain_def_points[1:], 1):
+ subdomain = mshr.Polygon(subdomain_vertices)
+ domain.set_subdomain(num, subdomain)
+ # count number of subdomains for further for loops
+ self._number_of_subdomains = num
# create global mesh. Submeshes still need to be extracted.
mesh = mshr.generate_mesh(domain, mesh_resolution)
+
if self.write2file['meshes_and_markers']:
filepath = self.output_dir+"global_mesh.pvd"
df.File(filepath) << mesh
filepath = self.output_dir+"global_mesh.xml.gz"
df.File(filepath) << mesh
+
self.mesh_subdomain.append(mesh)
# define domainmarker on global mesh and set marker values to domain
# number as assigned in the above for loop.
self.domain_marker = df.MeshFunction('size_t', mesh, mesh.topology().dim(), mesh.domains())
- # extract subdomain meshes
- for num in range(1, self._number_of_subdomains+1):
- # if MPI.size(mpi_comm_world()) == 1:
- self.mesh_subdomain.append(df.SubMesh(mesh, self.domain_marker, num))
- # else:
- # self.mesh_subdomain.append(create_submesh(mesh, self.domain_marker, num))
+ if self.interface_def_points is None:
+ self.domain_marker.set_all(0)
+ else:
+ # extract subdomain meshes
+ for num in range(1, self._number_of_subdomains+1):
+ # if MPI.size(mpi_comm_world()) == 1:
+ self.mesh_subdomain.append(df.SubMesh(mesh, self.domain_marker, num))
+ # else:
+ # self.mesh_subdomain.append(create_submesh(mesh, self.domain_marker, num))
@@ -884,21 +1040,25 @@ class LDDsimulation(object):
# define global interface marker. This is a mesh function on edges of the mesh.
self.interface_marker = df.MeshFunction('size_t', mesh, mesh.topology().dim()-1)
self.interface_marker.set_all(0)
- for glob_interface_num, vertices in enumerate(interface_def_points):
- self.interface.append(bi.Interface(vertices=vertices, #
- tol=self.tol,#mesh.hmin()/100,
- internal=True,#
- adjacent_subdomains = adjacent_subdomains[glob_interface_num],#
- isRichards = self.isRichards,#
- global_index = glob_interface_num,
- )
- )#
- # marker value gets internally set to global_index+1.
- # The reason is that we want to mark outer boundaries with the value 0.
- marker_value = self.interface[glob_interface_num].marker_value
- self.interface[glob_interface_num].mark(self.interface_marker, marker_value)
- self.interface[glob_interface_num].init_facet_dictionaries(self.interface_marker, marker_value, subdomain_index=0)
- # print(self.interface[glob_interface_num])
+ # if self.interface_def_points is None, i.e. we have no interface we do
+ # nothing, since self.init() defines self.interface = None
+ if self.interface_def_points is not None:
+ self.interface = []
+ for glob_interface_num, vertices in enumerate(interface_def_points):
+ self.interface.append(bi.Interface(vertices=vertices, #
+ tol=self.tol,#mesh.hmin()/100,
+ internal=True,#
+ adjacent_subdomains = adjacent_subdomains[glob_interface_num],#
+ isRichards = self.isRichards,#
+ global_index = glob_interface_num,
+ )
+ )#
+ # marker value gets internally set to global_index+1.
+ # The reason is that we want to mark outer boundaries with the value 0.
+ marker_value = self.interface[glob_interface_num].marker_value
+ self.interface[glob_interface_num].mark(self.interface_marker, marker_value)
+ self.interface[glob_interface_num].init_facet_dictionaries(self.interface_marker, marker_value, subdomain_index=0)
+ # print(self.interface[glob_interface_num])
if self.write2file['meshes_and_markers']:
filepath = self.output_dir+"interface_marker_global.pvd"
@@ -916,7 +1076,10 @@ class LDDsimulation(object):
# The list is_richards has the same length as the nuber of subdomains.
# Therefor it is used in the subdomain initialisation loop.
for subdom_num, isR in self.isRichards.items():
- interface_list = self._get_interfaces(subdom_num)
+ if self.interface_def_points is None:
+ interface_list = None
+ else:
+ interface_list = self._get_interfaces(subdom_num)
# print(f"has_interface list for subdomain {subdom_num}:", interface_list)
if self.densities is None:
subdom_densities = None
@@ -946,22 +1109,23 @@ class LDDsimulation(object):
)
# setup the linear solvers and set the solver parameters.
# df.LinearSolver(self.solver)
- if self.solver_type_is_Kryov:
+ if self.solver_type_is_Krylov:
self.subdomain[subdom_num].linear_solver = df.KrylovSolver(self.solver, self.preconditioner)
- # we use the previous iteration as an initial guess for the linear solver.
- solver_param = self.subdomain[subdom_num].linear_solver.parameters
- df.info(solver_param, True)
- solver_param['nonzero_initial_guess'] = self.solver_parameters['nonzero_initial_guess']
- solver_param['absolute_tolerance'] = self.solver_parameters['absolute_tolerance']
- solver_param['relative_tolerance'] = self.solver_parameters['relative_tolerance']
- solver_param['maximum_iterations'] = self.solver_parameters['maximum_iterations']
- solver_param['report'] = self.solver_parameters['report']
else:
self.subdomain[subdom_num].linear_solver = df.LUSolver()
+ print("solver method:")
+ for parameter, value in self.subdomain[subdom_num].linear_solver.parameters.items():
+ print(f"\'{parameter}\': {value}")
+
+ # assign solver parameters set at the beginning of the LDDsimulation class.
+ solver_param = self.subdomain[subdom_num].linear_solver.parameters
+ df.info(solver_param, False)
+ for parameter_name in self.solver_parameters.keys():
+ solver_param[parameter_name] = self.solver_parameters[parameter_name]
- # ## print out set solver parameters
- for parameter, value in self.subdomain[subdom_num].linear_solver.parameters.items():
- print(f"parameter: {parameter} = {value}")
+ # # ## print out set solver parameters
+ # for parameter, value in self.subdomain[subdom_num].linear_solver.parameters.items():
+ # print(f"parameter: {parameter} = {value}")
if self.write2file['meshes_and_markers']:
filepath = self.output_dir+f"interface_marker_subdomain{subdom_num}.pvd"
@@ -1030,12 +1194,30 @@ class LDDsimulation(object):
subdomain.pressure[phase].assign(pa0_interpolated)
subdomain.pressure_prev_iter[phase].assign(pa0_interpolated)
subdomain.pressure_prev_timestep[phase].assign(pa0_interpolated)
+ # write zeros to the abs_diff variables
+
+ absolute_difference = df.Function(subdomain.function_space["pressure"][phase])
+ # pressure_difference = pressure_exact[phase].vector()[:] - subdomain.pressure[phase].vector()[:]
+ # abs_diff_tmp = np.fabs(pressure_difference)
+ # absolute_difference.vector()[:] = abs_diff_tmp
+ if phase == "wetting":
+ data_set_name="abs_diff_w"
+ else:
+ data_set_name="abs_diff_nw"
+ self.write_function_to_xdmf(#
+ function=absolute_difference,
+ file=self.solution_file[num], #
+ time=self.starttime,#
+ data_set_name=data_set_name
+ )
# populate interface dictionaries with pressure values. The calc_gli_term
# of each subdomain reads from the interface.pressure_values_prev
# and from interface.pressure_values dictionary so both need to be populated.
- subdomain.write_pressure_to_interfaces()
- subdomain.write_pressure_to_interfaces(previous_iter = True)
+ if self.interface_def_points is not None:
+ subdomain.write_pressure_to_interfaces()
+ subdomain.write_pressure_to_interfaces(previous_iter = True)
+
subdomain.write_solution_to_xdmf(#
file = self.solution_file[num], #
time = self.starttime,#
@@ -1109,7 +1291,7 @@ class LDDsimulation(object):
interpolation_degree = self.interpolation_degree
if time is None:
- time = self.starttime
+ time = self.t
subdomain = self.subdomain[subdomain_index]
# in case time == self.starttime, the expression has already been
# assembled by self._init_DirichletBC_dictionary.
@@ -1128,6 +1310,65 @@ class LDDsimulation(object):
)
+ # def post_processing(self):
+ # """post processing of the simulation.
+ # calculate
+ # - pc_exact and pc_num
+ # - absolute differences to exact solution if present.
+ # - relative errors to exact solution if present
+ # """
+ # for subdom_ind, subdomain in self.subdomain.items():
+ # self._mesh = df.Mesh()
+ # solution_file = self.solution_file[subdom_ind]
+ # post_solution_file = self.postprocessed_solution_file[subdom_ind]
+ # pressure = dict()
+ # # the internal time series counter in df.XDMFFile starts at 0 and
+ # # refers to the first saved value. This corresponds to the initial
+ # # values. We then calculated self.number_of_timesteps many timesteps.
+ # t = self.starttime
+ # for internal_timestep in range(self.number_of_timesteps+1):
+ # for phase in subdomain.has_phases:
+ # pressure.update(
+ # {phase: df.Function(subdomain.function_space["pressure"][phase])}
+ # )
+ # phase_pressure_string = "pressure_"+"{}".format(phase)
+ # solution_file.read_checkpoint(
+ # u=pressure[phase],
+ # name=phase_pressure_string,
+ # counter=internal_timestep
+ # )
+ # pressure[phase].rename("pressure_"+"{}".format(phase), "pressure_"+"{}".format(phase))
+ # post_solution_file.write(pressure[phase], t)
+ # t += self.timestep_size
+ # print(f"read pressure of {phase}phase:")
+ # print(pressure[phase].vector()[:])
+ # solution_file.close()
+ # # # if we have an exact solution, write out the |u -uh|_L2 to see the
+ # # # absolute error.
+ # # if subdomain.pressure_exact is not None:
+ # # relative_L2_errornorm = dict()
+ # # for phase in subdomain.has_phases:
+ # # pa_exact = subdomain.pressure_exact[phase]
+ # # 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)
+ # # if norm_pa_exact > self.tol:
+ # # relative_L2_errornorm.update(
+ # # {phase: error_calculated/norm_pa_exact}
+ # # )
+ # # else:
+ # # relative_L2_errornorm.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 = relative_L2_errornorm,
+ # # )
+
+
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
@@ -1195,7 +1436,7 @@ class LDDsimulation(object):
if self.number_of_timesteps_to_analyse == 0:
self.analyse_timesteps = None
elif self.number_of_timesteps_to_analyse == 1:
- self.analyse_timesteps = [0]
+ self.analyse_timesteps = [1]
else:
if not self.number_of_timesteps%(self.number_of_timesteps_to_analyse-1) == 0:
# if number_of_timesteps%(number_of_timesteps_to_analyse-1) is not 0, we
@@ -1205,7 +1446,7 @@ class LDDsimulation(object):
new_analysing_interval = int(round(self.number_of_timesteps/(self.number_of_timesteps_to_analyse-1)))
self.number_of_timesteps = new_analysing_interval*(self.number_of_timesteps_to_analyse)
print(f"I'm resetting number_of_timesteps = {self.number_of_timesteps}\n")
- self.analyse_timesteps = np.linspace(0, self.number_of_timesteps, self.number_of_timesteps_to_analyse, dtype=int)
+ self.analyse_timesteps = np.linspace(1, self.number_of_timesteps, self.number_of_timesteps_to_analyse, dtype=int)
print(f"the following timesteps will be analysed: {self.analyse_timesteps}")
self.endtime = self.number_of_timesteps*self.timestep_size
diff --git a/LDDsimulation/boundary_and_interface.py b/LDDsimulation/boundary_and_interface.py
index b50d921..9204d77 100644
--- a/LDDsimulation/boundary_and_interface.py
+++ b/LDDsimulation/boundary_and_interface.py
@@ -333,7 +333,8 @@ class Interface(BoundaryPart):
interface_marker_value: int,
debug: bool = False) -> tp.Dict[int, tp.Dict[int, np.ndarray]]:
""" for a given function space function_space, calculate the dof indices
- and coordinates along the interface.
+ and coordinates along the interface. The pair {global_dof_index: dof_coordinate}
+ is saved as value in a dictionary which has global facet indices as keys.
# Parameters
function_space: #type df.functionSpace the function space of whitch
@@ -362,9 +363,21 @@ class Interface(BoundaryPart):
for mesh_entity in mesh_entity_iterator:
# this is the facet index relative to the Mesh
# that interface_marker is defined on.
- facet_index = mesh_entity.index()
+ interface_facet_index = mesh_entity.index()
+ for interface_facet in df.facets(mesh_entity):
+ # we only need the x,y coordinates. vertex.point()[:] returns
+ # the array and cannot be sliced. the result can be indexed, hence,
+ # vertex.point()[:][0:2]
+ # facet_midpoint = interface_facet.midpoint()[:][0:2]
+ # print(f"facet: {interface_facet}")
+ # print("facet coordinates: ")
+ facet_points = []
+ for vertex in df.vertices(interface_facet):
+ # print(f"vertex: {vertex} with coordinates {vertex.point()[:][0:2]}")
+ facet_points.append(vertex.point())
+
- dofs_and_coords.update({facet_index: dict()})
+ dofs_and_coords.update({interface_facet_index: dict()})
# because mesh_entity doesn't have a method to access coordinates we
# need to cast it into an acutal facet object.
for cell in df.cells(mesh_entity):
@@ -373,12 +386,15 @@ class Interface(BoundaryPart):
# returns an iterator object.
facet_cell = cell
facet_cell_index = facet_cell.index()
+ # cell_facets_entity = df.facets(cell)
+ # cell_coordinates = cell.get_vertex_coordinates()
# the cell is a triangle and for each dof numbered locally for that
- # cell, we have the coordinates in a numpy array.
+ # cell, we have the coordinates in a numpy array. The numbering is the
+ # numbering local to the cell.
facet_cell_dof_coordinates = element.tabulate_dof_coordinates(facet_cell)
# the cell dof map gives to each local (to the cell) numbering of the
- # dofs the index of the dof realitve to the global dof numbering of
+ # dofs the index of the dofs realitve to the global dof numbering of
# space, i.e. the index to access the dof with function.vector()[index].
facet_cell_dofmap = dofmap.cell_dofs(facet_cell_index)
@@ -387,9 +403,9 @@ class Interface(BoundaryPart):
# cell is a triangle and only has one facet that belongs to
# the interface. we want to store only dofs that lie on that particular
# facet
- if self._is_on_boundary_part(dof_coord):
+ if self._is_on_line_segment(facet_points[0], facet_points[1], dof_coord):
global_dof_ind = facet_cell_dofmap[local_dof_ind]
- dofs_and_coords[facet_index].update({global_dof_ind: dof_coord})
+ dofs_and_coords[interface_facet_index].update({global_dof_ind: dof_coord})
if debug:
@@ -482,10 +498,11 @@ class Interface(BoundaryPart):
for local_facet_index, local_facet_vertex in subdom_facet2coord.items():
for global_facet_index, global_facet_vertex in global_facet2coord.items():
# vertices are not necessarily in the same order
- vertices_are_equal = np.array_equal(local_facet_vertex[0], global_facet_vertex[0])
- vertices_are_equal += np.array_equal(local_facet_vertex[0], global_facet_vertex[1])
- vertices_are_equal += np.array_equal(local_facet_vertex[1], global_facet_vertex[0])
- vertices_are_equal += np.array_equal(local_facet_vertex[1], global_facet_vertex[1])
+ # np.array_equal(local_facet_vertex[0], global_facet_vertex[0])
+ vertices_are_equal = np.allclose(local_facet_vertex[0], global_facet_vertex[0], rtol=1e-10, atol=1e-14)
+ vertices_are_equal += np.allclose(local_facet_vertex[0], global_facet_vertex[1], rtol=1e-10, atol=1e-14)
+ vertices_are_equal += np.allclose(local_facet_vertex[1], global_facet_vertex[0], rtol=1e-10, atol=1e-14)
+ vertices_are_equal += np.allclose(local_facet_vertex[1], global_facet_vertex[1], rtol=1e-10, atol=1e-14)
both_vertices_are_equal = (vertices_are_equal == 2)
if both_vertices_are_equal:
self.local_to_global_facet_map[subdomain_index].update(
diff --git a/LDDsimulation/domainPatch.py b/LDDsimulation/domainPatch.py
index acdb82f..0b11785 100644
--- a/LDDsimulation/domainPatch.py
+++ b/LDDsimulation/domainPatch.py
@@ -209,25 +209,25 @@ class DomainPatch(df.SubDomain):
self._init_dof_maps()
self._init_markers()
self._init_measures()
- self._calc_interface_dof_indices_and_coordinates()
- self._init_interface_values()
- # self._calc_dof_indices_of_interfaces()
- # self._calc_interface_has_common_dof_indices()
- # self._calc_normal_vectors_norms_for_common_dofs()
+
if self.include_gravity:
self.gravity_term = dict()
self._calc_gravity_expressions()
else:
self.gravity_term = None
- # Dictionary of FEM function holding the same dofs gli_assembled above
- # but as a df.Function. This is needed for writing out the gli terms when
- # analyising the convergence of a timestep.
- self.gli_function = dict()
- for phase in self.has_phases:
- self.gli_function.update(
- {phase: df.Function(self.function_space["gli"][phase])}
- )
+ if self.has_interface is not None:
+ self._calc_interface_dof_indices_and_coordinates()
+ self._calc_corresponding_dof_indices()
+ self._init_interface_values()
+ # Dictionary of FEM function holding the same dofs gli_assembled above
+ # but as a df.Function. This is needed for writing out the gli terms when
+ # analyising the convergence of a timestep.
+ self.gli_function = dict()
+ for phase in self.has_phases:
+ self.gli_function.update(
+ {phase: df.Function(self.function_space["gli"][phase])}
+ )
# END constructor
@@ -238,16 +238,18 @@ class DomainPatch(df.SubDomain):
model = "two-phase"
if self.isRichards:
model = "Richards"
-
- interface_string = ["\nMy interfaces are:"]
- for glob_if_index in self.has_interface:
- interface = self.interface[glob_if_index]
- interface_string.append(interface.__str__())
- # common_int_dof_str = ".\n Dofs common with other interfaces:"\
- # +"{}".format([(ind, dof) for ind, dof in self._interface_has_common_dof_indices[glob_if_index].items()])
- # interface_string.append(common_int_dof_str)
- #
- interface_string = " ".join(interface_string)
+ if self.has_interface is not None:
+ interface_string = ["\nI have interfaces:"]
+ for glob_if_index in self.has_interface:
+ interface = self.interface[glob_if_index]
+ interface_string.append(interface.__str__())
+ # common_int_dof_str = ".\n Dofs common with other interfaces:"\
+ # +"{}".format([(ind, dof) for ind, dof in self._interface_has_common_dof_indices[glob_if_index].items()])
+ # interface_string.append(common_int_dof_str)
+ #
+ interface_string = " ".join(interface_string)
+ else:
+ interface_string ="\nI have no interfaces"
if self.outer_boundary is None:
outer_boundary_str = "no outer boundaries"
@@ -255,7 +257,6 @@ class DomainPatch(df.SubDomain):
outer_boundary_str = f"outer boundaries: {[ind for ind in self.outer_boundary_def_points.keys()]}.\n"
string = f"\nI'm subdomain{self.subdomain_index} and assume {model} model."\
- +f"\nI have the interfaces: {self.has_interface} and "\
+ outer_boundary_str\
+ interface_string
return string
@@ -273,22 +274,23 @@ class DomainPatch(df.SubDomain):
interface.pressure_values_prev
"""
subdomain = self.subdomain_index
- for ind in self.has_interface:
- interface = self.interface[ind]
- for phase in self.has_phases:
- p_dof_coordinates = self.interface_dof_indices_and_coordinates["pressure"][ind][phase]
- local2global_facet = interface.local_to_global_facet_map[subdomain]
- for subdomain_facet_index, facet_dof_coord_dict in p_dof_coordinates.items():
- global_facet_index = local2global_facet[subdomain_facet_index]
- # set dictionary to save the values to.
- if previous_iter:
- save_dict = interface.pressure_values_prev[subdomain][phase][global_facet_index]
- else:
- save_dict = interface.pressure_values[subdomain][phase][global_facet_index]
+ if self.has_interface is not None:
+ for ind in self.has_interface:
+ interface = self.interface[ind]
+ for phase in self.has_phases:
+ p_dof_coordinates = self.interface_dof_indices_and_coordinates["pressure"][ind][phase]
+ local2global_facet = interface.local_to_global_facet_map[subdomain]
+ for subdomain_facet_index, facet_dof_coord_dict in p_dof_coordinates.items():
+ global_facet_index = local2global_facet[subdomain_facet_index]
+ # set dictionary to save the values to.
+ if previous_iter:
+ save_dict = interface.pressure_values_prev[subdomain][phase][global_facet_index]
+ else:
+ save_dict = interface.pressure_values[subdomain][phase][global_facet_index]
- for dof_index, dof_coord in facet_dof_coord_dict.items():
- dof_coord_tupel = (dof_coord[0], dof_coord[1])
- save_dict.update({dof_coord_tupel: self.pressure[phase].vector()[dof_index]})
+ for dof_index, dof_coord in facet_dof_coord_dict.items():
+ dof_coord_tupel = (dof_coord[0], dof_coord[1])
+ save_dict.update({dof_coord_tupel: self.pressure[phase].vector()[dof_index]})
def governing_problem(self, phase: str) -> tp.Dict[str, fl.Form]:
@@ -341,17 +343,19 @@ class DomainPatch(df.SubDomain):
# and by interface_form terms if the neighbour of the current patch
# assumes Richards model
if phase == "wetting":
- # gather interface forms
- interface_forms = []
- gli_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")
- 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))
- gli = self.calc_gli_term(interface_index=interface, phase=phase)
- gli_forms.append((dt*gli*phi_a)*ds(marker_value))
+ if self.has_interface is not None:
+ # gather interface forms
+ interface_forms = []
+ gli_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")
+ 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))
+ gli = self.calc_gli_term(interface_index=interface, phase=phase)
+ gli_forms.append((dt*gli*phi_a)*ds(marker_value))
+
# form2 is different for wetting and nonwetting
form2 = dt/mu_a*df.dot(ka(S(pc_prev_iter))*df.grad(pa), df.grad(phi_a))*dx
rhs2 = -(porosity*(S(pc_prev_iter) - S(pc_prev_timestep))*phi_a)*dx
@@ -360,41 +364,50 @@ class DomainPatch(df.SubDomain):
rhs_gravity = dt/mu_a*df.dot(ka(S(pc_prev_iter))*df.grad(z_a), df.grad(phi_a))*dx
elif phase == "nonwetting":
- # gather interface forms
- interface_forms = []
- gli_forms = []
- for interface in self.has_interface:
- marker_value = self.interface[interface].marker_value
- if self.interface[interface].neighbour_isRichards[self.subdomain_index]:
- # if the neighbour of our subdomain (with index self.subdomain_index)
- # assumes a Richards model, we assume constant normalised atmospheric pressure
- # and no interface term is practically appearing.
- # interface_forms += (df.Constant(0)*phi_a)*ds(interface)
- # pass
- interface_forms.append((0*pa*phi_a)*ds(marker_value))
- gli_forms.append((0*phi_a)*ds(marker_value))
- else:
- # interface_forms += (dt*Lambda_a*pa*phi_a)*ds(interface)
- interface_forms.append((dt*Lambda_a*pa*phi_a)*ds(marker_value))
- gli = self.calc_gli_term(interface_index=interface, phase=phase)
- gli_forms.append((dt*gli*phi_a)*ds(marker_value))
+ if self.has_interface is not None:
+ # gather interface forms
+ interface_forms = []
+ gli_forms = []
+ for interface in self.has_interface:
+ marker_value = self.interface[interface].marker_value
+ if self.interface[interface].neighbour_isRichards[self.subdomain_index]:
+ # if the neighbour of our subdomain (with index self.subdomain_index)
+ # assumes a Richards model, we assume constant normalised atmospheric pressure
+ # and no interface term is practically appearing.
+ # interface_forms += (df.Constant(0)*phi_a)*ds(interface)
+ # pass
+ interface_forms.append((0*pa*phi_a)*ds(marker_value))
+ gli_forms.append((0*phi_a)*ds(marker_value))
+ else:
+ # interface_forms += (dt*Lambda_a*pa*phi_a)*ds(interface)
+ interface_forms.append((dt*Lambda_a*pa*phi_a)*ds(marker_value))
+ gli = self.calc_gli_term(interface_index=interface, phase=phase)
+ gli_forms.append((dt*gli*phi_a)*ds(marker_value))
+
form2 = dt/mu_a*df.dot(ka(1 - S(pc_prev_iter))*df.grad(pa), df.grad(phi_a))*dx
# the non wetting phase has a + before instead of a - before the bracket
rhs2 = (porosity*(S(pc_prev_iter) - S(pc_prev_timestep))*phi_a)*dx
if self.include_gravity:
- # z_a included the minus sign already, see self._calc_gravity_expressions
+ # z_a included the minus sign, see self._calc_gravity_expressions
rhs_gravity = dt/mu_a*df.dot(ka(1 - S(pc_prev_iter))*df.grad(z_a), df.grad(phi_a))*dx
else:
raise RuntimeError('missmatch for input parameter phase.',
'Expected either phase == wetting or phase == nonwetting ')
return None
- form = form1 + form2 + sum(interface_forms)
- if self.include_gravity:
- # z_a included the minus sign already, see self._calc_gravity_expressions
- rhs = rhs1 + rhs2 + dt*source_a*phi_a*dx - sum(gli_forms) - rhs_gravity
+ if self.has_interface is not None:
+ form = form1 + form2 + sum(interface_forms)
+ if self.include_gravity:
+ # z_a included the minus sign already, see self._calc_gravity_expressions
+ rhs = rhs1 + rhs2 + dt*source_a*phi_a*dx - sum(gli_forms) + rhs_gravity
+ else:
+ rhs = rhs1 + rhs2 + dt*source_a*phi_a*dx - sum(gli_forms)
else:
- rhs = rhs1 + rhs2 + dt*source_a*phi_a*dx - sum(gli_forms)
-
+ form = form1 + form2
+ if self.include_gravity:
+ # z_a included the minus sign already, see self._calc_gravity_expressions
+ rhs = rhs1 + rhs2 + dt*source_a*phi_a*dx + rhs_gravity
+ else:
+ rhs = rhs1 + rhs2 + dt*source_a*phi_a*dx
form_and_rhs = {#
'form': form,#
'rhs' : rhs
@@ -455,52 +468,43 @@ class DomainPatch(df.SubDomain):
# to calculate a gl0 term.
# TODO: add intrinsic permeabilty!!!!
if self.include_gravity:
- # z_a included the minus sign already, see self._calc_gravity_expressions
- flux = -1/mu_a*ka(S(pc_prev))*df.grad(pa_prev + z_a)
+ # z_a included the minus sign, see self._calc_gravity_expressions
+ flux = -1/mu_a*ka(S(pc_prev))*df.grad(pa_prev - z_a)
else:
flux = -1/mu_a*ka(S(pc_prev))*df.grad(pa_prev)
else:
# phase == 'nonwetting'
# we need anothre flux in this case
if self.include_gravity:
- flux = -1/mu_a*ka(1-S(pc_prev))*df.grad(pa_prev + z_a)
+ flux = -1/mu_a*ka(1-S(pc_prev))*df.grad(pa_prev - z_a)
else:
flux = -1/mu_a*ka(1-S(pc_prev))*df.grad(pa_prev)
flux_function = df.project(flux, self.function_space["flux"][phase])
flux_x, flux_y = flux_function.split()
gl0 = df.Function(self.function_space["gli"][phase])
-
+ # # get dictionaries of global dof indices and coordinates along
+ # # the interface. These dictionaries have the facet indices of
+ # # facets belonging to the interface as keys (with respect to
+ # # the submesh numbering) and the dictionary
+ # # containing pairs {global_dof_index: dof_coordinate} for all
+ # # dofs along the facet as values.
gli_dof_coordinates = self.interface_dof_indices_and_coordinates["gli"][interf_ind][phase]
- p_dof_coordinates = self.interface_dof_indices_and_coordinates["pressure"][interf_ind][phase]
- flux_dof_coordinates = self.interface_dof_indices_and_coordinates["flux"][interf_ind][phase]
-
+ # p_dof_coordinates = self.interface_dof_indices_and_coordinates["pressure"][interf_ind][phase]
+ # flux_dof_coordinates = self.interface_dof_indices_and_coordinates["flux"][interf_ind][phase]
+ # # the attributes local and global for facet numbering mean
+ # # the following: local facet index is the index of the facet
+ # # with respect to the numbering of the submesh belonging to
+ # # the subdomain. global facet index refers to the facet numbering
+ # # of the mesh of the whole computational domain.
local2global_facet = interface.local_to_global_facet_map[subdomain]
+ corresponding_dof_index = self.interface_corresponding_dof_index[interf_ind][phase]
for local_facet_ind, normal in interface.outer_normals[subdomain].items():
gli_dofs_along_facet = gli_dof_coordinates[local_facet_ind]
- p_dofs_along_facet = p_dof_coordinates[local_facet_ind]
- flux_x_dofs_along_facet = flux_dof_coordinates["x"][local_facet_ind]
- flux_y_dofs_along_facet = flux_dof_coordinates["y"][local_facet_ind]
-
- # for each gli dof with index gli_dof_ind, we need the
- # dof indices for the flux and the pressure corresponding
- # to the coordinates of the dof, in order to calculate gl0
- # dof wise.
for gli_dof_index, gli_dof_coord in gli_dofs_along_facet.items():
- flux_x_dof_index = -1
- for flux_x_dof_ind, dof_coord in flux_x_dofs_along_facet.items():
- if np.array_equal(dof_coord, gli_dof_coord):
- flux_x_dof_index = flux_x_dof_ind
-
- flux_y_dof_index = -1
- for flux_y_dof_ind, dof_coord in flux_y_dofs_along_facet.items():
- if np.array_equal(dof_coord, gli_dof_coord):
- flux_y_dof_index = flux_y_dof_ind
-
- p_dof_index = -1
- for p_dof_ind, dof_coord in p_dofs_along_facet.items():
- if np.array_equal(dof_coord, gli_dof_coord):
- p_dof_index = p_dof_ind
+ flux_x_dof_index = corresponding_dof_index[local_facet_ind][gli_dof_index]["flux_x"]
+ flux_y_dof_index = corresponding_dof_index[local_facet_ind][gli_dof_index]["flux_y"]
+ p_dof_index = corresponding_dof_index[local_facet_ind][gli_dof_index]["pressure"]
gl0.vector()[gli_dof_index] = flux_x.vector()[flux_x_dof_index]*normal[0] \
+ flux_y.vector()[flux_y_dof_index]*normal[1] \
@@ -587,11 +591,11 @@ class DomainPatch(df.SubDomain):
"and iteration = " + "{} on subdomain".format(iteration)+ "{}".format(subdomain)\
+ "You suck! Revisite DomainPatch.calc_gli_term() and LDDsimulation.LDDsolver ASAP.")
- dt = self.timestep_size
+ # dt = self.timestep_size
Lambda = self.lambda_param[phase]
gli_interface_dofs = self.interface_dof_indices_and_coordinates["gli"][interface_index][phase]
- p_interface_dofs = self.interface_dof_indices_and_coordinates["pressure"][interface_index][phase]
+ # p_interface_dofs = self.interface_dof_indices_and_coordinates["pressure"][interface_index][phase]
local2global_facet = interface.local_to_global_facet_map[subdomain]
for local_facet_index, gli_facet_dofs2coordinates in gli_interface_dofs.items():
global_facet_index = local2global_facet[local_facet_index]
@@ -599,20 +603,42 @@ class DomainPatch(df.SubDomain):
gli_prev = interface.gli_term_prev[neighbour][phase][global_facet_index]
p_prev = pressure_prev[global_facet_index]
gli_save_dict = gli_save_dictionary[global_facet_index]
+ # print(f"\non subdomain{self.subdomain_index}")
+ # print(f"on facet with glob_index: {global_facet_index} we have gli_dofs")
+ #
for gli_dof_index, gli_dof_coord in gli_facet_dofs2coordinates.items():
# find the right previous pressure dof of the interface
# dictionary to use for the calculation.
- p_previous_dof = -9999999.9999
+ # print(f"coordinates of gli_dof with index {gli_dof_index} I try to match is: {gli_dof_coord}" )
+ p_previous_dof = None
for p_dof_coord_tupel, p_prev_dof in p_prev.items():
p_coord = np.array([p_dof_coord_tupel[0], p_dof_coord_tupel[1]])
- if np.array_equal(gli_dof_coord, p_coord):
+ # print(f"coordinates of p_dof I test this against: {p_coord}")
+ a_close_b = np.allclose(gli_dof_coord, p_coord, rtol=1e-10, atol=1e-14)
+ b_close_a = np.allclose(p_coord, gli_dof_coord, rtol=1e-10, atol=1e-14)
+ # print(f"type(gli_dof_coord) = {type(gli_dof_coord)}")
+ # print(f"np.allclose({gli_dof_coord}, {p_coord}, rtol=1e-10, atol=1e-14): {a_close_b }")
+ # print(f"np.allclose({p_coord}, {gli_dof_coord}, rtol=1e-10, atol=1e-14): {b_close_a}")
+ if a_close_b and b_close_a:
p_previous_dof = p_prev_dof
+ print("\n")
+ if p_previous_dof is None:
+ raise RuntimeError(f"didn't find p_previous_dof corresponding to dict key ({gli_dof_coord})",
+ "something is wrong")
# same with the gli_previous_dof
- gli_previous_dof = -9999999.9999
+ gli_previous_dof = None
for gli_prev_coord_tupel, gli_prev_dof in gli_prev.items():
gli_prev_coord = np.array([gli_prev_coord_tupel[0], gli_prev_coord_tupel[1]])
- if np.array_equal(gli_dof_coord, gli_prev_coord):
+ # due to the warning in the documentation of np.allclose
+ # that np.allclose is not symetric, we test if both
+ # np.allclose(a,b) and np.allclose(b,a) are true.
+ a_close_b = np.allclose(gli_dof_coord, gli_prev_coord, rtol=1e-10, atol=1e-14)
+ b_close_a = np.allclose(gli_prev_coord, gli_dof_coord, rtol=1e-10, atol=1e-14)
+ if a_close_b and b_close_a:
gli_previous_dof = gli_prev_dof
+ if gli_previous_dof is None:
+ raise RuntimeError(f"didn't find gli_previous_dof corresponding to dict key ({gli_dof_coord})",
+ "something is wrong")
gli_value = -2*Lambda*p_previous_dof - gli_previous_dof
gli.vector()[gli_dof_index] = gli_value
@@ -641,8 +667,11 @@ class DomainPatch(df.SubDomain):
of the subdomain.
"""
subdomain = self.subdomain_index
- for index in self.has_interface:
- self.interface[index].current_iteration[subdomain] = 0
+ if self.has_interface is not None:
+ for index in self.has_interface:
+ self.interface[index].current_iteration[subdomain] = 0
+ else:
+ pass
#### PRIVATE METHODS
def _init_function_space(self) -> None:
@@ -681,13 +710,13 @@ class DomainPatch(df.SubDomain):
""" calculate dof maps and the vertex to parent map.
This method sets the class Variables
- self.parent_mesh_index
self.dof2vertex
self.vertex2dof
"""
mesh_data = self.mesh.data()
- # local to parent vertex index map
- self.parent_mesh_index = mesh_data.array('parent_vertex_indices',0)
+ # # local to parent vertex index map
+ # if self.has_interface is not None:
+ # 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}")
self.dof2vertex = dict()
self.vertex2dof = dict()
@@ -720,6 +749,7 @@ class DomainPatch(df.SubDomain):
{"y": self.function_space["flux"][phase].sub(1).dofmap()}
)
+
def _init_markers(self) -> None:
""" define boundary markers
@@ -732,17 +762,18 @@ class DomainPatch(df.SubDomain):
self.outer_boundary_marker.set_all(0)
self.interface_marker = df.MeshFunction('size_t', self.mesh, self.mesh.topology().dim()-1)
self.interface_marker.set_all(0)
- for glob_index in self.has_interface:
- # the marker value is set to the same as the global marker value,
- # stored in self.interfaces[glob_index].marker_value.
- marker_value = self.interface[glob_index].marker_value
- # each interface gets marked with the global interface index.
- self.interface[glob_index].mark(self.interface_marker, marker_value)
- # init facet_to_vertex_coordinates_dictionray for subdoaain mesh
- self.interface[glob_index].init_facet_dictionaries(
- interface_marker=self.interface_marker,
- interface_marker_value=marker_value,
- subdomain_index=self.subdomain_index)
+ if self.has_interface is not None:
+ for glob_index in self.has_interface:
+ # the marker value is set to the same as the global marker value,
+ # stored in self.interfaces[glob_index].marker_value.
+ marker_value = self.interface[glob_index].marker_value
+ # each interface gets marked with the global interface index.
+ self.interface[glob_index].mark(self.interface_marker, marker_value)
+ # init facet_to_vertex_coordinates_dictionray for subdoaain mesh
+ self.interface[glob_index].init_facet_dictionaries(
+ interface_marker=self.interface_marker,
+ interface_marker_value=marker_value,
+ subdomain_index=self.subdomain_index)
# create outerBoundary objects and mark outer boundaries with 1
# in case there is no outer boundary, self.outer_boundary_def_points is
# None.
@@ -766,6 +797,7 @@ class DomainPatch(df.SubDomain):
self.outer_boundary_marker, boundary_marker_value
)
+
def _init_measures(self):
""" define measures for the governing form
@@ -814,6 +846,91 @@ class DomainPatch(df.SubDomain):
)
+ def _calc_corresponding_dof_indices(self):
+ """ calculate dictionary which for each interface and each phase holds
+ for each facet index, the dof indices of the pressures and flux components
+ corresponding to a given dof index of the gli function.
+ This gets used by self.calc_gl0_term
+ """
+ self.interface_corresponding_dof_index = dict()
+ subdomain = self.subdomain_index
+ for interf_ind in self.has_interface:
+ interface = self.interface[interf_ind]
+ self.interface_corresponding_dof_index.update(
+ {interf_ind: dict()}
+ )
+ for phase in self.has_phases:
+ self.interface_corresponding_dof_index[interf_ind].update(
+ {phase: dict()}
+ )
+ # get dictionaries of global dof indices and coordinates along
+ # the interface. These dictionaries have the facet indices of
+ # facets belonging to the interface as keys (with respect to
+ # the submesh numbering) and the dictionary
+ # containing pairs {global_dof_index: dof_coordinate} for all
+ # dofs along the facet as values.
+ gli_dof_coordinates = self.interface_dof_indices_and_coordinates["gli"][interf_ind][phase]
+ p_dof_coordinates = self.interface_dof_indices_and_coordinates["pressure"][interf_ind][phase]
+ flux_dof_coordinates = self.interface_dof_indices_and_coordinates["flux"][interf_ind][phase]
+ # the attributes local and global for facet numbering mean
+ # the following: local facet index is the index of the facet
+ # with respect to the numbering of the submesh belonging to
+ # the subdomain. global facet index refers to the facet numbering
+ # of the mesh of the whole computational domain.
+ # local2global_facet = interface.local_to_global_facet_map[subdomain]
+ for local_facet_ind in interface.outer_normals[subdomain].keys():
+ self.interface_corresponding_dof_index[interf_ind][phase].update(
+ {local_facet_ind: dict()}
+ )
+ gli_dofs_along_facet = gli_dof_coordinates[local_facet_ind]
+ p_dofs_along_facet = p_dof_coordinates[local_facet_ind]
+ flux_x_dofs_along_facet = flux_dof_coordinates["x"][local_facet_ind]
+ flux_y_dofs_along_facet = flux_dof_coordinates["y"][local_facet_ind]
+
+ # for each gli dof with index gli_dof_ind, we need the
+ # dof indices for the flux and the pressure corresponding
+ # to the coordinates of the dof, in order to calculate gl0
+ # dof wise.
+ for gli_dof_index, gli_dof_coord in gli_dofs_along_facet.items():
+ self.interface_corresponding_dof_index[interf_ind][phase][local_facet_ind].update(
+ {gli_dof_index: {"flux_x": None,
+ "flux_x": None,
+ "pressure": None}}
+ )
+ for flux_x_dof_ind, dof_coord in flux_x_dofs_along_facet.items():
+ a_close_b = np.allclose(gli_dof_coord, dof_coord, rtol=1e-10, atol=1e-14)
+ b_close_a = np.allclose(dof_coord, gli_dof_coord, rtol=1e-10, atol=1e-14)
+ if a_close_b and b_close_a:
+ self.interface_corresponding_dof_index[interf_ind][phase][local_facet_ind][gli_dof_index].update(
+ {"flux_x": flux_x_dof_ind}
+ )
+ if self.interface_corresponding_dof_index[interf_ind][phase][local_facet_ind][gli_dof_index]["flux_x"] is None:
+ raise RuntimeError(f"didn't find flux_x_dof_ind corresponding to dict key ({gli_dof_coord})",
+ "something is wrong")
+
+ for flux_y_dof_ind, dof_coord in flux_y_dofs_along_facet.items():
+ a_close_b = np.allclose(gli_dof_coord, dof_coord, rtol=1e-10, atol=1e-14)
+ b_close_a = np.allclose(dof_coord, gli_dof_coord, rtol=1e-10, atol=1e-14)
+ if a_close_b and b_close_a:
+ self.interface_corresponding_dof_index[interf_ind][phase][local_facet_ind][gli_dof_index].update(
+ {"flux_y": flux_y_dof_ind}
+ )
+ if self.interface_corresponding_dof_index[interf_ind][phase][local_facet_ind][gli_dof_index]["flux_y"] is None:
+ raise RuntimeError(f"didn't find flux_y_dof_ind corresponding to dict key ({gli_dof_coord})",
+ "something is wrong")
+
+ for p_dof_ind, dof_coord in p_dofs_along_facet.items():
+ a_close_b = np.allclose(gli_dof_coord, dof_coord, rtol=1e-10, atol=1e-14)
+ b_close_a = np.allclose(dof_coord, gli_dof_coord, rtol=1e-10, atol=1e-14)
+ if a_close_b and b_close_a:
+ self.interface_corresponding_dof_index[interf_ind][phase][local_facet_ind][gli_dof_index].update(
+ {"pressure": p_dof_ind}
+ )
+ if self.interface_corresponding_dof_index[interf_ind][phase][local_facet_ind][gli_dof_index]["pressure"] is None:
+ raise RuntimeError(f"didn't find p_dof_ind corresponding to dict key ({gli_dof_coord})",
+ "something is wrong")
+
+
def _init_interface_values(self):
""" allocate the dictionaries for each interface used to store the
interface values for communication.
@@ -860,17 +977,24 @@ class DomainPatch(df.SubDomain):
write_iter_for_fixed_time bool flag to toggle wether pairs of
(solution, iteration) for fixed t should be
written instead of (solution, t)
- subsequent_errors df.Funtion FEM Function passed by the solver
- containing self.pressure - self.pressure_prev_iter
"""
t = time
- for phase in self.has_phases:
- if not write_iter_for_fixed_time:
+ if not write_iter_for_fixed_time:
+ for phase in self.has_phases:
self.pressure[phase].rename("pressure_"+"{}".format(phase), "pressure_"+"{}".format(phase))
file.write(self.pressure[phase], t)
+ capillary_pressure = df.Function(self.function_space["pressure"]["wetting"])
+ if self.isRichards:
+ capillary_pressure.assign(-self.pressure["wetting"])
else:
+ pc_temp = self.pressure["nonwetting"].vector()[:] - self.pressure["wetting"].vector()[:]
+ capillary_pressure.vector()[:] = pc_temp
+ capillary_pressure.rename("pc_num", "pc_num")
+ file.write(capillary_pressure, t)
+ else:
+ for phase in self.has_phases:
i = self.iteration_number
self.pressure[phase].rename("pressure_"+"{}".format(phase), #
"pressure(t="+"{}".format(t)+")_"+"{}".format(phase)
@@ -883,9 +1007,5 @@ class DomainPatch(df.SubDomain):
"pressure(t="+"{}".format(t)+")_"+"{}".format(phase)
)
file.write(error, i)
- # self.gli_function[phase].vector().set_local(self.gli_assembled[phase])
- # self.gli_function[phase].rename("gli_function_"+"{}".format(phase),#
- # "all_gli_terms(t="+"{}".format(t)+")_"+"{}".format(phase))
- # file.write(self.gli_function[phase], i)
# END OF CLASS DomainPatch
diff --git a/LDDsimulation/helpers.py b/LDDsimulation/helpers.py
index 95c6168..765630c 100644
--- a/LDDsimulation/helpers.py
+++ b/LDDsimulation/helpers.py
@@ -8,6 +8,9 @@ import sys
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
+import typing as tp
+import sympy as sym
+
# from termcolor import colored
# from colormap import Colormap
#import fenicstools as ft
@@ -18,174 +21,112 @@ def print_once(*args,**kwargs):
print(*args,**kwargs)
sys.stdout.flush()
-# def fit_IC(mesh,u):
-# """
-# Fit given data u on mesh to certain functions.
-# Useful in the case of a single droplet.
-# Returns the fitted data as class for IC.
-# """
-#
-# # get mesh coordiantes
-# xx = mesh.coordinates()[:,0]
-# yy = mesh.coordinates()[:,1]
-# nx = np.size(xx)
-# ny = np.size(yy)
-#
-# # get function values
-# rho_data = np.zeros(nx)
-# phase_field_data = np.zeros(nx)
-# mu_data = np.zeros(nx)
-#
-# for i in range(0,nx):
-# rho_data[i] = u.sub(0)(xx[i],yy[i])
-# mu_data[i] = u.sub(4)(xx[i],yy[i])
-# phase_field_data[i] = u.sub(2)(xx[i],yy[i])
-#
-# # define functions to fit
-# def rho(x,a,b,cx,cy,dx,dy):
-# val = a*np.tanh(b*np.sqrt((x[0]-cx)**2+(x[1]-cy)**2)+dx) + dy
-# return val
-#
-# def phase_field(x,a,b,cx,cy,dx,dy):
-# val = a*np.tanh(b*np.sqrt((x[0]-cx)**2+(x[1]-cy)**2)+dx) + dy
-# return val
-#
-# def mu(x,a,b,cx,cy,dx,dy):
-# val = a*np.tanh(b*np.sqrt((x[0]-cx)**2+(x[1]-cy)**2)+dx) + dy
-# return val
-#
-# xx1 = np.concatenate((xx[:,np.newaxis],yy[:,np.newaxis]),axis= 1)
-# # fit data to function
-# popt_rho = curve_fit(rho, xx1.T , rho_data, p0=[1./np.pi,-10,.5,.3,1.,.5])
-# [a,b,cx,cy,dx,dy] = popt_rho[0]
-# # print(a,b,cx,cy,dx,dy)
-#
-# popt_phase_field = curve_fit(phase_field, xx1.T , phase_field_data, p0=[1./np.pi,-10,.5,.3,1.,.5])
-# [a1,b1,cx1,cy1,dx1,dy1] = popt_phase_field[0]
-# # print(a1,b1,cx1,cy1,dx1,dy1)
-#
-# popt_mu = curve_fit(mu, xx1.T , mu_data, p0=[1./np.pi,-10,.5,.3,1.,.5])
-# [a2,b2,cx2,cy2,dx2,dy2] = popt_mu[0]
-# # print(a2,b2,cx2,cy2,dx2,dy2)
-#
-# # plot result
-# # res = 500;
-# # [X,Y] = np.meshgrid(np.linspace(0.,1.,res),np.linspace(0.,1.,res))
-# # Z = np.zeros((res,res))
-# # for i in range(0,res):
-# # for j in range(0,res):
-# # Z[i,j] = phase_field([X[i,j],Y[i,j]],a1,b1,cx1,cy1,d1)
-# # plt.imshow(Z,cmap = "viridis", origin ="lower")
-# # plt.colorbar()
-# # plt.show()
-#
-# class fitted_IC(Expression):
-# def eval(self, values, x):
-# # rho
-# values[0] = rho(x,a,b,cx,cy,dx,dy)
-#
-# # velocity
-# for i in range(2):
-# values[i + 1] = 0.
-# # c
-# values[2 + 1] = phase_field(x,a1,b1,cx1,cy1,dx1,dy1)
-# # tau
-# values[2 + 2] = 0.0
-# # mu
-# values[2 + 3] = mu(x,a2,b2,cx2,cy2,dx2,dy2)
-# # sigma
-# for i in range(2):
-# values[2 + 4 + i] = 0.
-#
-# def value_shape(self):
-# return (4 + 2*2,)
-#
-# return fitted_IC
-#
-#
-# def tensor_jump(v, n):
-# """
-# DG jump operator for vector valued functions
-# """
-# return df.outer(v("+"), n("+")) + df.outer(v("-"), n("-"))
-#
-# def print_constants(self,file=None):
-# """
-# Print the names and values of all df.Constant attributes.
-# """
-# for key in self.__dict__.keys():
-# attribute = self.__dict__[key]
-# if type(attribute) is type(df.Constant(0.)):
-# print_once(attribute.name() + " = " + str(attribute.values()),file=file)
-# for key in self.__dict__["_EOS"].__dict__:
-# attribute = self.__dict__["_EOS"].__dict__[key]
-# if type(attribute) is type(df.Constant(0.)):
-# print_once(attribute.name() + " = " + str(attribute.values()),file=file)
-#
-# def get_interface_thickness(solution,ME,xmin=-.5,xmax=1.5,ymin=0.,ymax=1.,ny=1000):
-# """
-# Compute interface thickness from simulation result.
-# The interface thickness is here defined as the region where 0.1 < c < 0.9.
-# Parameter: solution, function space ME
-# Meshrelated: xmin,xmax,ymin,ymax
-# Number of generated samples: ny
-# """
-# xmid = .5*(xmin+xmax)
-# yy = np.linspace(ymin,ymax,ny)
-# xx = xmid*np.ones_like(yy)
-# # evaluate phasefield along following points
-# zz = np.column_stack((xx,yy))
-# p = ft.Probes(zz.flatten(), ME.sub(2).collapse())
-# p(solution.sub(2))
-# # values collected on root process
-# values= p.array(root=0)
-# if df.MPI.rank(df.MPI.comm_world) == 0:
-# try:
-# # get interface_start
-# interface_start = np.min(np.where(abs(values-0.1)<5.e-3))
-# # get interface_end
-# interface_end = np.min(np.where(abs(values-0.9)<5.e-3))
-# # compute interface_thickness
-# interface_thickness = abs(yy[interface_end]-yy[interface_start])
-# return interface_thickness
-# except:
-# print_once(colored("Phasefield values out of range.","red"))
-# return -1
-#
-#
-# def plot_solution(u):
-# """ Plot solution for poster or the like..."""
-# # Define "CD"-Colormap
-# red = [0.12549019607843137, 0.3843137254901961, 0.6431372549019608, 0.9019607843137255]
-# blue = [0.6823529411764706, 0.788235294117647, 0.8941176470588236, 1.0]
-# green = [0.33725490196078434, 0.5254901960784314, 0.7137254901960784, 0.9019607843137255]
-# c = Colormap()
-# mycmap = c.cmap( {'red':red[::-1], 'green':green[::-1], 'blue':blue[::-1]})
-# matplotlib.cm.register_cmap(name='mycmap', cmap=mycmap)
-# # plt.rcParams['image.cmap'] = 'mycmap'
-# plt.rcParams['image.cmap'] = 'RdYlBu'
-# matplotlib.rcParams.update({'font.size': 17})
-# matplotlib.rc('font',**{'family':'sans-serif','sans-serif':['Helvetica']})
-# matplotlib.rc('text', usetex=True)
-#
-# cnt = plot(u,mode="contourf",levels=256,antialiased = True,linestyles=None)
-# # fixes to remove white contour lines
-# cnt = plot(u,mode="contourf",levels=256,antialiased = True,linestyles=None)
-# cnt = plot(u,mode="contourf",levels=256,antialiased = True,linestyles=None)
-#
-# for c in cnt.collections:
-# c.set_edgecolor("face")
-#
-# plt.xlim(0.45, 0.8)
-# plt.ylim(0.2, 0.55)
-# plt.axis("off")
-# plt.annotate("$+$ phase",color='w',xy=(0.5, .45))
-# plt.annotate("$-$ phase",xy=(0.7, .3))
-# plt.annotate("$\Omega$",xy=(0.75, .5),fontsize=25)
-# plt.annotate("$\gamma$",xy=(0.647, .358),fontsize=25)
-# plt.annotate("",xy=(0.62,0.375),xytext=(0.66,0.33),arrowprops=dict(arrowstyle="<->",
-# connectionstyle="arc3",facecolor='black',linewidth=1))
-# #plt.show()
-# plt.tight_layout()
-# plt.savefig('diffuse_simulation.pdf',bbox_inches='tight')
-# return 0
+
+def generate_exact_solution_expressions(
+ symbols: tp.Dict[int, int],
+ isRichards: tp.Dict[int, bool],
+ symbolic_pressure: tp.Dict[int, tp.Dict[str, str]],
+ symbolic_capillary_pressure: tp.Dict[int, str],
+ viscosity: tp.Dict[int, tp.List[float]],#
+ porosity: tp.Dict[int, tp.Dict[str, float]],#
+ relative_permeability: tp.Dict[int, tp.Dict[str, tp.Callable[...,None]] ],#
+ relative_permeability_prime: tp.Dict[int, tp.Dict[str, tp.Callable[...,None]] ],
+ saturation_pressure_relationship: tp.Dict[int, tp.Callable[...,None]] = None,#
+ saturation_pressure_relationship_prime: tp.Dict[int, tp.Callable[...,None]] = None,#
+ symbolic_S_pc_relationship: tp.Dict[int, tp.Callable[...,None]] = None,#
+ symbolic_S_pc_relationship_prime: tp.Dict[int, tp.Callable[...,None]] = None,#
+ densities: tp.Dict[int, tp.Dict[str, float]] = None,#
+ gravity_acceleration: float = 9.81,
+ include_gravity: bool = False,
+ ) -> tp.Dict[str, tp.Dict[int, tp.Dict[str, str]] ]:
+ """ Generate exact solution, source and initial condition code from symbolic expressions"""
+
+ output = {
+ 'source': dict(),
+ 'initial_condition': dict(),
+ 'exact_solution': dict()
+ }
+ x = symbols["x"]
+ y = symbols["y"]
+ t = symbols["t"]
+ # turn above symbolic code into exact solution for dolphin and
+ # construct the rhs that matches the above exact solution.
+ dtS = dict()
+ div_flux = dict()
+ if saturation_pressure_relationship is not None:
+ S_pc_sym = saturation_pressure_relationship
+ S_pc_sym_prime = saturation_pressure_relationship_prime
+ for subdomain, isR in isRichards.items():
+ dtS.update({subdomain: dict()})
+ div_flux.update({subdomain: dict()})
+ output['source'].update({subdomain: dict()})
+ output['exact_solution'].update({subdomain: dict()})
+ output['initial_condition'].update({subdomain: dict()})
+ if isR:
+ subdomain_has_phases = ["wetting"]
+ else:
+ subdomain_has_phases = ["wetting", "nonwetting"]
+
+ # conditional for S_pc_prime
+ pc = symbolic_capillary_pressure[subdomain]
+ dtpc = sym.diff(pc, t, 1)
+ dxpc = sym.diff(pc, x, 1)
+ dypc = sym.diff(pc, y, 1)
+ if saturation_pressure_relationship is not None:
+ S = sym.Piecewise((S_pc_sym[subdomain](pc), pc > 0), (1, True))
+ dS = sym.Piecewise((S_pc_sym_prime[subdomain](pc), pc > 0), (0, True))
+ else:
+ S = symbolic_S_pc_relationship[subdomain]
+ dS = symbolic_S_pc_relationship_prime[subdomain]
+ for phase in subdomain_has_phases:
+ # Turn above symbolic expression for exact solution into c code
+ output['exact_solution'][subdomain].update(
+ {phase: sym.printing.ccode(symbolic_pressure[subdomain][phase])}
+ )
+ # save the c code for initial conditions
+ output['initial_condition'][subdomain].update(
+ {phase: sym.printing.ccode(symbolic_pressure[subdomain][phase].subs(t, 0))}
+ )
+ if phase == "nonwetting":
+ dtS[subdomain].update(
+ {phase: -porosity[subdomain]*dS*dtpc}
+ )
+ else:
+ dtS[subdomain].update(
+ {phase: porosity[subdomain]*dS*dtpc}
+ )
+ pa = symbolic_pressure[subdomain][phase]
+ dxpa = sym.diff(pa, x, 1)
+ dxdxpa = sym.diff(pa, x, 2)
+ dypa = sym.diff(pa, y, 1)
+ dydypa = sym.diff(pa, y, 2)
+ mu = viscosity[subdomain][phase]
+ ka = relative_permeability[subdomain][phase]
+ dka = relative_permeability_prime[subdomain][phase]
+ rho = densities[subdomain][phase]
+ g = gravity_acceleration
+
+ if phase == "nonwetting":
+ # x part of div(flux) for nonwetting
+ dxdxflux = -1/mu*dka(1-S)*dS*dxpc*dxpa + 1/mu*dxdxpa*ka(1-S)
+ # y part of div(flux) for nonwetting
+ if include_gravity:
+ dydyflux = -1/mu*dka(1-S)*dS*dypc*(dypa + rho*g) + 1/mu*dydypa*ka(1-S)
+ else:
+ dydyflux = -1/mu*dka(1-S)*dS*dypc*dypa + 1/mu*dydypa*ka(1-S)
+ else:
+ # x part of div(flux) for wetting
+ dxdxflux = 1/mu*dka(S)*dS*dxpc*dxpa + 1/mu*dxdxpa*ka(S)
+ # y part of div(flux) for wetting
+ if include_gravity:
+ dydyflux = 1/mu*dka(S)*dS*dypc*(dypa + rho*g) + 1/mu*dydypa*ka(S)
+ else:
+ dydyflux = 1/mu*dka(S)*dS*dypc*(dypa) + 1/mu*dydypa*ka(S)
+
+ div_flux[subdomain].update({phase: dxdxflux + dydyflux})
+ contructed_rhs = dtS[subdomain][phase] - div_flux[subdomain][phase]
+ output['source'][subdomain].update(
+ {phase: sym.printing.ccode(contructed_rhs)}
+ )
+
+ return output
diff --git a/LDDsimulation/solutionFile.py b/LDDsimulation/solutionFile.py
index 7ae0a95..e71450f 100644
--- a/LDDsimulation/solutionFile.py
+++ b/LDDsimulation/solutionFile.py
@@ -9,5 +9,5 @@ class SolutionFile(df.XDMFFile):
df.XDMFFile.__init__(self, mpicomm, filepath)
self.parameters["functions_share_mesh"] = True
- self.parameters["flush_output"] = False
+ self.parameters["flush_output"] = True
self.path = filepath # Mimic the file path attribute from a `file` returned by `open`
--
GitLab