diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index a9d2afc4c958b4e9c1eddc3c42d65a57035d52ea..01ca7f6912fae264874a971a0bf23914222f4427 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -9,6 +9,8 @@ import typing as tp
import mshr
import domainPatch as dp
import helpers as hlp
+import pandas as pd
+
class LDDsimulation(object):
""" Main Class of LDD simulation
@@ -19,7 +21,9 @@ class LDDsimulation(object):
def __init__(self, dimension: int = 2,#
tol: float = None):
-
+ hlp.print_once("\n ###############################################\n",
+ "############# Begin Simulation ################\n",
+ "###############################################\n")
hlp.print_once("Init LDD simulation...")
# parameter as set in NSAC_simulation by Lukas Ostrowski
hlp.print_once("Set internal fenics parameter...")
@@ -34,14 +38,18 @@ class LDDsimulation(object):
#Parameters
# df.parameters["allow_extrapolation"] = True
# df.parameters["refinement_algorithm"] = "plaza_with_parent_facets"
- df.parameters["form_compiler"]["quadrature_degree"] = 10
+ df.parameters["form_compiler"]["quadrature_degree"] = 12
+ # interpolation degree, for source terms, intitial and boundary conditions.
+ self.interpolation_degree = 8
# # To be able to run DG in parallel
# df.parameters["ghost_mode"] = "shared_facet"
+ # df.parameters["ghost_mode"] = "none"
# 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.INFO)
# df.set_log_level(df.LogLevel.INFO)
@@ -53,7 +61,7 @@ class LDDsimulation(object):
# TRACE = 13, // what's happening (in detail)
# DBG = 10 // sundry
- hlp.print_once("Set default parameter...")
+ hlp.print_once("Set default parameter...\n")
self.restart_file = None
# # write checkpoints
# self.write_checkpoints = True
@@ -61,6 +69,10 @@ class LDDsimulation(object):
# self.project_to_cg = False
# # no sponge zones
# self.sponge = False
+ # set up mpi for internal parallelisation.
+ self._mpi_rank = df.MPI.rank(df.MPI.comm_world)
+ self._mpi_communicator = df.MPI.comm_world
+
## Public variables.
# directory to write files output to.
@@ -81,13 +93,15 @@ class LDDsimulation(object):
## Private variables
# maximal number of L-iterations that the LDD solver uses.
- self._max_iter_num = 100
+ self._max_iter_num = 300
# TODO rewrite this with regard to the mesh sizes
self.calc_tol = self.tol
# The number of subdomains are counted by self.init_meshes_and_markers()
self._number_of_subdomains = 0
# variable to check if self.set_parameters() has been called
self._parameters_set = False
+ # flag holds true if self.initialise() has been executed.
+ self._initialised = False
# dictionary of objects of class DomainPatch initialised by self._init_subdomains()
self.subdomain = dict()
# dictionary to hold the input Expressions for the external boundary. These are
@@ -111,6 +125,9 @@ class LDDsimulation(object):
# These are interpolated and stored to each respective subdomain by
# self._init_initial_values()
self.initial_conditions = None
+ # dictionary of filenames for writing the solution to xdmf files for each
+ # subdomain. This gets populated by self._init_solution_files()
+ self.solution_filename = None
########## SOLVER DEFINITION AND PARAMETERS ########
# print("\nLinear Algebra Backends:")
# df.list_linear_algebra_backends()
@@ -120,8 +137,19 @@ class LDDsimulation(object):
# df.list_krylov_solver_preconditioners()
### Define the linear solver to be used.
- self.solver = 'bicgstab' # biconjugate gradient stabilized method
- self.preconditioner = 'hypre_amg' # incomplete LU factorization
+ self.solver = 'bicgstab' #'gmres'#'bicgstab' # biconjugate gradient stabilized method
+ self.preconditioner = '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-12,
+ 'relative_tolerance': 1E-10,
+ 'maximum_iterations': 1000,
+ 'report': True
+ }
+
def set_parameters(self,#
output_dir: str,#
@@ -138,6 +166,7 @@ class LDDsimulation(object):
lambda_param: tp.Dict[int, tp.Dict[str, float]],#
relative_permeability: tp.Dict[int, tp.Dict[str, tp.Callable[...,None]] ],#
saturation: tp.Dict[int, tp.Callable[...,None]],#
+ time_interval: tp.List[float],#
timestep_size: tp.Dict[int, float],#
sources: tp.Dict[int, tp.Dict[str, str]],#
initial_conditions: tp.Dict[int, tp.Dict[str, str]],#
@@ -160,6 +189,10 @@ class LDDsimulation(object):
self.lambda_param = lambda_param
self.relative_permeability = relative_permeability
self.saturation = saturation
+ # simulation start time and time variable
+ self.t = time_interval[0]
+ # end_time
+ self.T = time_interval[1]
self.timestep_size = timestep_size
self.sources = sources
self.initial_conditions = initial_conditions
@@ -175,18 +208,66 @@ class LDDsimulation(object):
"""
if not self._parameters_set:
raise RuntimeError('Parameters note set!\nYou forgott to run self.set_parameters(**kwds)')
+ print("initialise meshes and marker functions ...\n")
self._init_meshes_and_markers()
+ print("initialise interfaces ...\n")
self._init_interfaces()
+ print("initialise subdomains ...\n")
self._init_subdomains()
+ # initial values get initialised here to be able to call the solver at different
+ # timesteps without having to call self.run().
self._init_initial_values()
+ self._initialised = True
+
+ def run(self, analyse_solver_at_times: tp.List[float] = None):
+ """ run the simulation
+
+ This method implements the time stepping, calls the acutal solvers and writes
+ solutions to files. It is the main funtion of the simulation class. If you
+ want to understand this code, start reading here.
- def LDDsolver(self, time: float, debug: bool = False):
+ PARAMETERS
+ analyse_solver_at_times #type tp.List[float] List of timesteps at which
+ to produce output to
+ analyse the solver at
+ given timestep. This
+ tells the function at
+ which timesteps to
+ toggle the analyse_timestep
+ flag of self.LDDsolver.
+
+ """
+ timer = df.Timer("Simulation time")
+ # check if simulation has been initialised. If the user calls self.initialise()
+ # from his/her script, then the time for initialisation is not counted in
+ # the simulation time.
+ if not self._initialised:
+ print("Simulation has not been initialised by the user. \n",
+ "Initialising simulation now ... ")
+ self.initialise()
+
+
+ #
+
+ def LDDsolver(self, time: float = None, #
+ debug: bool = False, #
+ analyse_timestep: bool = False):
""" calulate the solution on all patches for the next timestep
The heart of the simulation. The function LDDsolver implements the LDD
iteration, that is solves the problem given the initial data at timestep
time.
"""
+ # in case no time is given, use the time of the class.
+ if time is None:
+ time = self.t
+
+ if analyse_timestep:
+ # dictionary holding filename strings for each subdomain to store
+ # iterations in, used to analyse the behaviour of the iterates within
+ # the current time step.
+ solution_over_iteration_within_timestep = dict()
+ subsequent_errors_over_iteration = dict()
# dictionary holding bool variables for each subdomain indicating wether
# minimal error has been achieved, i.e. the calculation can be considered
# done or not.
@@ -194,20 +275,18 @@ class LDDsimulation(object):
# before the iteration gets started all iteration numbers must be nulled
# and the subdomain_calc_isdone dictionary must be set to False
for ind, subdomain in self.subdomain.items():
- #create solver object for subdomains
- if time < self.calc_tol:
- subdomain.linear_solver = df.KrylovSolver(self.solver, self.preconditioner)
- # we use the previous iteration as an initial guess for the linear solver.
- solver_param = subdomain.linear_solver.parameters
- solver_param['nonzero_initial_guess'] = True
- # solver_param['absolute_tolerance'] = 1E-12
- # solver_param['relative_tolerance'] = 1E-9
- solver_param['maximum_iterations'] = 1000
- solver_param['report'] = True
- # ## print out set solver parameters
- # for parameter, value in self.linear_solver.parameters.items():
- # print(f"parameter: {parameter} = {value}")
- # df.info(solver_param, True)
+ # create xdmf files for writing out iterations if analyse_timestep is
+ # given.
+ if analyse_timestep:
+ filename = self.output_dir+"subdomain"+"{}_".format(ind)+ "iterations_at_time"+ "{}".format(time) +".xdmf"
+ 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()} )
+ #
subdomain_calc_isdone.update({ind: False})
# reset all interface[has_interface].current_iteration[ind] = 0, for
@@ -215,10 +294,11 @@ class LDDsimulation(object):
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(interpolation_degree = 3,
- time = time,
+ self._eval_sources(time = time,
subdomain_index = ind)
- self.update_DirichletBC_dictionary(time = time, subdomain_index = ind)
+ self.update_DirichletBC_dictionary(time = time, #
+ subdomain_index = ind,
+ )
# the following only needs to be done when time is not equal 0 since
# our initialisation functions already took care of setting what follows
# for the initial iteration step at time = 0.
@@ -247,32 +327,67 @@ class LDDsimulation(object):
# set subdomain iteration to current iteration
subdomain.iteration_number = iteration
if debug:
- print(f"Entering iteration {iteration}")
+ print(f"On subdomain{sd_index}: Entering iteration {iteration}\n")
+ print(f"subdomain{sd_index}.isRichards = {subdomain.isRichards}.")
+ print(f"subdomain{sd_index}.has_phases = {subdomain.has_phases}.")
# solve the problem on subdomain
self.Lsolver_step(subdomain_index = sd_index,#
debug = debug,
)
- subsequent_iter_error = dict()
- for phase in subdomain.has_phases:
- pa = subdomain.pressure[phase]
- pa_prev_i = subdomain.pressure_prev_iter[phase]
- subsequent_iter_error.update(
- {phase: df.errornorm(pa, pa_prev_i, 'L2')}
- )
-
- total_subsequent_iter_error = max(subsequent_iter_error.values())
- # check if error criterion has been met
- if total_subsequent_iter_error < subdomain.mesh.hmin()/60:
- subdomain_calc_isdone[sd_index] = True
+ if iteration > 0:
+ subsequent_iter_error = dict()
+ for phase in subdomain.has_phases:
+ error = df.Function(subdomain.function_space[phase])
+ error.assign(subdomain.pressure[phase] - subdomain.pressure_prev_iter[phase])
+ error_calculated = df.norm(error, 'L2')
+ subsequent_iter_error.update(
+ {phase: error_calculated}
+ )
+ print(f"the subsequent error in iteration {iteration} is {subsequent_iter_error[phase]}")
+
+ if analyse_timestep:
+ #
+ subsequent_errors_over_iteration[sd_index].update(
+ { phase: error })
+ subsequent_error_filename = self.output_dir+"subdomain"+"{}_".format(sd_index)+\
+ "subsequent_iteration_error_for_phase_"+"{}".format(phase)+"_at_time"+ "{}".format(time) +".csv"
+ self.write_subsequent_errors(
+ filename = subsequent_error_filename, #
+ subdomain_index = sd_index,
+ subsequent_errors = subsequent_iter_error,
+ time = time)
+
+
+ # both phases should be done before we mark the subdomain as
+ # finished.
+ total_subsequent_iter_error = max(subsequent_iter_error.values())
+ # check if error criterion has been met
+ error_criterion = subdomain.mesh.hmin()**2*subdomain.timestep_size
+ if debug:
+ print(f" total_subsequent_error in iter {iteration} = {total_subsequent_iter_error }\n")
+ print(f"the grid size on subdomain{sd_index} is h={subdomain.mesh.hmin()}")
+ print(f"the error criterion is tol = {error_criterion}")
+ if total_subsequent_iter_error < error_criterion:
+ if debug:
+ print(f"subdomain{sd_index} reachd error tol = {error_criterion} in iteration {iteration}.")
+ subdomain_calc_isdone[sd_index] = True
# prepare next iteration
# write the newly calculated solution to the inteface dictionaries
# for communication
subdomain.write_pressure_to_interfaces()
for phase in subdomain.has_phases:
+ if analyse_timestep:
+ subdomain.write_solution_to_xdmf(#
+ file = solution_over_iteration_within_timestep[sd_index], #
+ time = time,#
+ write_iter_for_fixed_time = True,
+ )#subsequent_errors = subsequent_errors_over_iteration[sd_index])
+
subdomain.pressure_prev_iter[phase].assign(
subdomain.pressure[phase]
)
+
else:
# one subdomain is done. Check if all subdomains are done to
# stop the while loop if needed.
@@ -340,7 +455,9 @@ class LDDsimulation(object):
# apply outer Dirichlet boundary conditions if present.
if subdomain.outer_boundary is not None:
- self.outerBC[subdomain_index][phase].apply(form_assembled, rhs_assembled)
+ # apply boundary dirichlet conditions for all outer boundaries.
+ for index, boundary in enumerate(subdomain.outer_boundary):
+ self.outerBC[subdomain_index][index][phase].apply(form_assembled, rhs_assembled)
if debug and subdomain.mesh.num_cells() < 36:
print("\nSystem after applying outer boundary conditions:")
@@ -364,6 +481,76 @@ class LDDsimulation(object):
print("\npressure after solver:\n", subdomain.pressure[phase].vector().get_local())
+ def _init_solution_files(self):
+ """ set up solution files for saving the solution of the LDD scheme for
+ each subdomain.
+ """
+ # set up solution file names for each subdomain.
+ self.solution_filename = dict()
+ for subdom_ind, subdomain in self.subdomain.items():
+ tau = subdomain.timestep_size
+ L = subodmain.L
+ h = subdomain.mesh.hmin()
+ if subdomain.isRichards:
+ Lam = subdomain.lambda_param['wetting']
+ filename = self.output_dir+"solution_subdomain" + "{}_".format(subdom_ind)#
+ +"dt" + "{}".format(tau)#
+ # only show three digits of the mesh size.
+ +"_hmin" + "{number:.{digits}f}".format(number=h, digits=3)#
+ +"_lambda" + "{}".format(Lam)
+ + "_Lw" + "{}".format(L["wetting"])#
+ +".xdmf"
+ else:
+ Lam_w = subdomain.lambda_param['wetting']
+ Lam_nw = subdomain.lambda_param['nonwetting']
+ filename = self.output_dir+"solution_subdomain" + "{}_".format(subdom_ind)#
+ +"dt" + "{}".format(tau)#
+ # only show three digits of the mesh size.
+ +"_hmin" + "{number:.{digits}f}".format(number=h, digits=3)#
+ +"_lambda_w" + "{}".format(Lam_w)#
+ +"_lambda_nw" + "{}".format(Lam_nw)#
+ +"_Lw" + "{}".format(L["wetting"])#
+ +"_Lnw" + "{}".format(L["nonwetting"])
+ +".xdmf"
+
+ self.solution_filename.update( #
+ {ind: SolutionFile(self._mpi_communicator, filename)}
+ )
+
+
+ def write_subsequent_errors(self,#
+ filename: str, #
+ subdomain_index: int,#
+ subsequent_errors,#: tp.Dict[str: float],#
+ time: float = None) -> None:
+ """ write subsequent errors to csv file for plotting with pgfplots"""
+ # in case no time is given, use the time of the class.
+ subdomain = self.subdomain[subdomain_index]
+ if time is None:
+ time = self.t
+
+ subdomain = self.subdomain[subdomain_index]
+ iteration = subdomain.iteration_number
+
+ if subdomain.isRichards:
+ self.errors = pd.DataFrame({ "iteration": iteration,
+ #"th. initial mass": df.assemble(rho*df.dx), \
+ #"energy": self.total_energy(), \
+ "wetting": subsequent_errors["wetting"]}, index=[iteration])
+
+ else:
+ self.errors = pd.DataFrame({"iteration":iteration, #\
+ #"th. initial mass": df.assemble(rho*df.dx), \
+ #"energy": self.total_energy(), \
+ "wetting": subsequent_errors["wetting"],
+ "nonwetting": subsequent_errors["nonwetting"]})
+ # if self._rank == 0:
+ # check if file exists
+ if os.path.isfile(filename) == True:
+ with open(filename, 'a') as f:
+ self.errors.to_csv(f, header=False, sep='\t', encoding='utf-8', index=False)
+ else:
+ self.errors.to_csv(filename, sep='\t', encoding='utf-8', index=False)
## Private methods
def _init_meshes_and_markers(self, subdomain_def_points: tp.List[tp.List[df.Point]] = None,#
@@ -429,7 +616,12 @@ class LDDsimulation(object):
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.write2file['meshes_and_markers']:
filepath = self.output_dir+"domain_marker.pvd"
@@ -489,20 +681,25 @@ class LDDsimulation(object):
internal=True,#
adjacent_subdomains = adjacent_subdomains[num],#
isRichards = self.isRichards,#
- global_index = num)
+ global_index = num,
+ )
)#
- # the marking number of each interface corresponds to the mathematical
- # numbering of the interfaces (starting with 1 not 0 for the first interface).
+ # marker value gets internally set to global_index+1.
# The reason is that we want to mark outer boundaries with the value 0.
- self.interface[num].mark(self.interface_marker, num+1)
- self.interface[num].init_interface_values(self.interface_marker, num+1)
+ self.interface[num].mark(self.interface_marker, self.interface[num].marker_value)
+ print("Global interface vertices:")
+ self.interface[num].init_interface_values(self.interface_marker, self.interface[num].marker_value)
if self.write2file['meshes_and_markers']:
filepath = self.output_dir+"interface_marker_global.pvd"
df.File(filepath) << self.interface_marker
def _init_subdomains(self) -> None:
- """ generate list of opjects of class DomainPatch.
+ """ generate dictionary of opjects of class DomainPatch.
+
+ This method generates a dictionary of opjects of class DomainPatch, containing
+ everything that is needed for the LDDsolver. Additionally a linear solver
+ object is created.
"""
# The list is_richards has the same length as the nuber of subdomains.
@@ -528,6 +725,20 @@ class LDDsimulation(object):
tol = self.tol)
}
)
+ # setup the linear solvers and set the solver parameters.
+ 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
+ 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']
+ # ## print out set solver parameters
+ # for parameter, value in self.linear_solver.parameters.items():
+ # print(f"parameter: {parameter} = {value}")
+ # df.info(solver_param, True)
+
@@ -555,12 +766,15 @@ class LDDsimulation(object):
return interface_of_subdomain
- def _init_initial_values(self, interpolation_degree: int = 2):
+ def _init_initial_values(self, interpolation_degree: int = None):
""" set initial values
This method populates the subdomain dictionaries as well as the interface
dictionaries with initial values needed for the LDDsolver iterations.
"""
+ if interpolation_degree is None:
+ interpolation_degree = self.interpolation_degree
+
for num, subdomain in self.subdomain.items():
mesh = subdomain.mesh
V = subdomain.function_space
@@ -585,10 +799,10 @@ class LDDsimulation(object):
subdomain.pressure_prev_timestep[phase].assign(pa0_interpolated)
# populate interface dictionaries with pressure values.
- subdomain.write_pressure_to_interfaces(initial_values = True)
+ subdomain.write_pressure_to_interfaces()
def update_DirichletBC_dictionary(self, subdomain_index: int,#
- interpolation_degree: int = 2, #
+ interpolation_degree: int = None, #
time: float = 0,#
):
""" update time of the dirichlet boundary object for subdomain with
@@ -597,6 +811,9 @@ class LDDsimulation(object):
In case we don't have a case with an exact solution, we should have 0
expressions in self.dirichletBC_Expressions, so nothing should happen.
"""
+ if interpolation_degree is None:
+ interpolation_degree = self.interpolation_degree
+
num = subdomain_index
subdomain = self.subdomain[num]
# we only need to set Dirichlet boundary conditions if we have an outer
@@ -604,36 +821,58 @@ class LDDsimulation(object):
if subdomain.outer_boundary is not None:
mesh = subdomain.mesh
V = subdomain.function_space
- # Here the dictionary has to be created first because t = 0.
- if np.abs(time) < self.calc_tol:
- self.outerBC.update({num: dict()})
- self.dirichletBC_dfExpression.update({num: dict()})
-
- pDC = self.dirichletBC_Expressions[num]
- boundary_marker = subdomain.outer_boundary_marker
- for phase in subdomain.has_phases:
- # note that the default interpolation degree is 2
- if np.abs(time) < self.tol :
- # time = 0 so the Dolfin Expression has to be created first.
- pDCa = df.Expression(pDC[phase], domain = mesh, degree = interpolation_degree, t=time)
- self.dirichletBC_dfExpression[num].update({phase: pDCa})
- else:
- pDCa = self.dirichletBC_dfExpression[num][phase].t = time
-
- self.outerBC[num].update(
- {phase: df.DirichletBC(V[phase], pDCa, boundary_marker, 1)}
- )
+ # loop over all outer boundary parts. This is a bit of a brainfuck:
+ # subdomain.outer_boundary gives a dictionary of the outer boundaries of the subdomain.
+ # Since a domain patch can have several disjoint outer boundary parts, the expressions
+ # need to get an enumaration index which starts at 0. So subdomain.outer_boundary[j] is
+ # the dictionary of outer dirichlet conditions of the subdomain and boundary part j for
+ # the wetting and the nonwetting phase, so subdomain.outer_boundary[j]['wetting'] and
+ # subdomain.outer_boundary[j]['nonwetting'] return
+ # the actual expression needed for the dirichlet condition for both phases if present.
+ for index, boundary in enumerate(subdomain.outer_boundary):
+ if np.abs(time) < self.tol:
+ # Here the dictionary has to be created first because t = 0.
+ # create a dictionary of dictionaries for each subdomain
+ self.outerBC.update({num: dict()})
+ self.dirichletBC_dfExpression.update({num: dict()})
+ # create a dictionary of dictionaries for each outer boundary
+ # part of the subdomain.
+ self.outerBC[num].update({index: dict()})
+ self.dirichletBC_dfExpression[num].update({index: dict()})
+
+ pDC = self.dirichletBC_Expressions[num][index]
+ boundary_marker = subdomain.outer_boundary_marker
+ # as with the interfaces, since numbering of the outer boundary
+ # parts of each subdomain starts at 0 and subdomain.outer_boundary_marker
+ # gets set to 0 on all facets of the mesh, we need to up the value for
+ # an actual outer boundary part by 1.
+ boundary_marker_value = index+1
+ for phase in subdomain.has_phases:
+ # note that the default interpolation degree is 2
+ if np.abs(time) < self.tol :
+ # time = 0 so the Dolfin Expression has to be created first.
+ pDCa = df.Expression(pDC[phase], domain = mesh, degree = interpolation_degree, t=time)
+ self.dirichletBC_dfExpression[num][index].update({phase: pDCa})
+ else:
+ pDCa = self.dirichletBC_dfExpression[num][index][phase].t = time
+
+ self.outerBC[num][index].update(
+ {phase: df.DirichletBC(V[phase], pDCa, boundary_marker, boundary_marker_value)}
+ )
else:
# case subdomain.outer_boundary is None
print("subdomain is an inner subdomain and has no outer boundary.\n",
"no Dirichlet boundary has been set.")
def _eval_sources(self, subdomain_index: int, #
- interpolation_degree: int = 2, #
+ interpolation_degree: int = None, #
time: float = 0,#
):
""" evaluate time dependent source terms or initialise them if time == 0
"""
+ if interpolation_degree is None:
+ interpolation_degree = self.interpolation_degree
+
subdomain = self.subdomain[subdomain_index]
for phase in subdomain.has_phases:
if np.abs(time) < self.calc_tol:
@@ -652,3 +891,17 @@ class LDDsimulation(object):
else:
# note that the default interpolation degree is 2
subdomain.source[phase].t = time
+
+
+# https://github.com/geo-fluid-dynamics/phaseflow-fenics
+# adapted from Lucas Ostrowski
+class SolutionFile(df.XDMFFile):
+ """
+ This class extends `df.XDMFFile` with some minor changes for convenience.
+ """
+ def __init__(self, mpicomm, filepath):
+ df.XDMFFile.__init__(self, mpicomm, filepath)
+
+ self.parameters["functions_share_mesh"] = True
+ self.parameters["flush_output"] = True
+ self.path = filepath # Mimic the file path attribute from a `file` returned by `open`
diff --git a/LDDsimulation/domainPatch.py b/LDDsimulation/domainPatch.py
index 1ba9132dfec044c7a4b6e68a7f897b87419e9cb0..6bb78942124064b50cdb358b65018f08d9dfbf41 100644
--- a/LDDsimulation/domainPatch.py
+++ b/LDDsimulation/domainPatch.py
@@ -144,7 +144,14 @@ class DomainPatch(df.SubDomain):
self.saturation = saturation
# timestep size, tau in the paper
self.timestep_size = timestep_size
+
### Class variables set by methods
+ # sometimes we need to loop over the phases and the length of the loop is
+ # determined by the model we are assuming
+ if self.isRichards:
+ self.has_phases =['wetting']
+ else:
+ self.has_phases =['wetting', 'nonwetting']
# dictionary holding the dof indices corresponding to an interface of
# given interface. see self._calc_dof_indices_of_interfaces()
self._dof_indices_of_interface = dict()
@@ -169,13 +176,6 @@ class DomainPatch(df.SubDomain):
self.pressure_prev_timestep = dict()
# test function(s) for the form set by _init_function_space
self.testfunction = None
- # Dictionary of Arrays (one for 'wetting' and one for 'nonwetting')
- # corresponding to the dofs of an FEM function which is initiated
- # by calc_gli_term() holding the dofs of the assembled dt*gli*v*ds terms as
- # a sparse vector (here, dt = self.timestep_size). This vector is read by the solver and added to the rhs
- # which is given by self.govering_problem().
- #
- self.gli_assembled = dict()
# 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().
@@ -184,24 +184,32 @@ class DomainPatch(df.SubDomain):
# first time LDDsimulation.LDDsolver is run.
self.linear_solver = None
- # sometimes we need to loop over the phases and the length of the loop is
- # determined by the model we are assuming
- if self.isRichards:
- self.has_phases =['wetting']
- else:
- self.has_phases =['wetting', 'nonwetting']
-
self._init_function_space()
- self._init_dof_and_vertex_maps()
+ self._init_dof_maps()
self._init_markers()
self._init_measures()
self._calc_dof_indices_of_interfaces()
self._calc_interface_has_common_dof_indices()
+ # Dictionary of Arrays (one for 'wetting' and one for 'nonwetting')
+ # corresponding to the dofs of an FEM function which is initiated
+ # by calc_gli_term() holding the dofs of the assembled dt*gli*v*ds terms as
+ # a sparse vector (here, dt = self.timestep_size). This vector is read by the solver and added to the rhs
+ # which is given by self.govering_problem().
+ #
+ self.gli_assembled = dict()
+ # 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[phase])}
+ )
# END constructor
#### PUBLIC METHODS
- def write_pressure_to_interfaces(self, initial_values: bool = False):
+ def write_pressure_to_interfaces(self):
""" save the interface values of self.pressure to all neighbouring interfaces
This method is used by the LDDsimulation.LDDsolver method.
@@ -212,12 +220,12 @@ class DomainPatch(df.SubDomain):
This is needed for the first ever iteration
at time = 0.
"""
- if initial_values:
- # this assumes, that self.pressure_prev_iter has been set by
- # LDDsimulation._init_initial_values().
- from_func = self.pressure_prev_iter
- else:
- from_func = self.pressure
+ # if initial_values:
+ # # this assumes, that self.pressure_prev_iter has been set by
+ # # LDDsimulation._init_initial_values().
+ # from_func = self.pressure_prev_iter
+ # else:
+ from_func = self.pressure
for ind in self.has_interface:
for phase in self.has_phases:
@@ -280,7 +288,8 @@ class DomainPatch(df.SubDomain):
interface_forms = []
# the marker values of the interfacese are global interfaces index + 1
for interface in self.has_interface:
- marker_value = interface + 1
+ 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))
# form2 is different for wetting and nonwetting
@@ -290,14 +299,14 @@ class DomainPatch(df.SubDomain):
# gather interface forms
interface_forms = []
for interface in self.has_interface:
- marker_value = interface + 1
+ 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((df.Constant(0)*phi_a)*ds(marker_value))
+ # pass
+ interface_forms.append((df.Constant(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))
@@ -317,7 +326,7 @@ class DomainPatch(df.SubDomain):
return form_and_rhs
- def calc_gli_term(self) -> None: # tp.Dict[str, fl.Form]
+ def calc_gli_term(self, debug: bool = False) -> None: # tp.Dict[str, fl.Form]
"""calculate the gl terms for the iteration number of the subdomain
calculate the gl terms for the iteration number of the subdomain and save them
@@ -344,12 +353,11 @@ class DomainPatch(df.SubDomain):
)
for ind in self.has_interface:
- marker_value = ind + 1
interface = self.interface[ind]
# neighbour index
neighbour = interface.neighbour[subdomain]
neighbour_iter_num = interface.current_iteration[neighbour]
- ds = self.ds(marker_value)
+ 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
@@ -357,7 +365,12 @@ class DomainPatch(df.SubDomain):
# by 1/2 to to get an average of the gli terms for dofs belonging to
# two different interfaces.
dofs_in_common_with_other_interfaces = self._interface_has_common_dof_indices[ind]
-
+ if debug:
+ print(f"On subdomain{subdomain}, we have:\n",
+ f"Interface{ind}, Neighbour index = {neighbour}\n",
+ f"dofs on interface:\n{interface_dofs['wetting']}\n",
+ f"dof2global_vertex_map:\n{self.parent_mesh_index[self.dof2vertex['wetting'][interface_dofs['wetting']]]}\n"
+ f"and dofs common with other interfaces:\n{dofs_in_common_with_other_interfaces['wetting']}")
if iteration == 0:
n = df.FacetNormal(self.mesh)
p_prev_iter = self.pressure_prev_timestep
@@ -368,9 +381,17 @@ class DomainPatch(df.SubDomain):
if not neighbour_iter_num == 0:
# copy the gl term from the neighbour first to use it in the
# next iteration
+ if debug:
+ 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()}
)
+
+ # for glob_dof_ind, dof in interface.gli_term[neighbour]['wetting'].items():
+ # interface.gli_term_prev[subdomain]['wetting'][glob_dof_ind] = dof
+ 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'])
# in case we have two-phase and the neighbour is two-phase, two,
# we need to Additionally save the nonwetting values.
if not interface.neighbour_isRichards[subdomain] and not self.isRichards:
@@ -408,6 +429,9 @@ class DomainPatch(df.SubDomain):
# there are no dofs that lie on another interface and we can
# skip the following step.
if dofs_in_common_with_other_interfaces[phase].size != 0:
+ if debug:
+ print(f"the following dofs (phase = {phase}) of interface{ind} are in common with \n",
+ f" other interfaces: {dofs_in_common_with_other_interfaces[phase]}")
for common_dof in dofs_in_common_with_other_interfaces[phase]:
# from the standpoint of the subdomain we are currently on,
# each dof can belong maximally to two interfaces. On these
@@ -429,8 +453,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
@@ -466,10 +490,17 @@ class DomainPatch(df.SubDomain):
# gli_prev_term of the current subdomain. This is because
# we assume, that the neighbouring patch has done a previous
# iteration and we need the value to calculate the gli term with
+ 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()}
)
+ 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'])
+
if not interface.neighbour_isRichards[subdomain] and not self.isRichards:
interface.gli_term_prev[subdomain].update(#
{'nonwetting': interface.gli_term[neighbour]['nonwetting'].copy()}
@@ -483,6 +514,10 @@ class DomainPatch(df.SubDomain):
# the gli_assembled array of the subdomain
# to be able to sum the gli terms of all interfaces together,
# we need a dummy vector to hold the dofs read from the interface.
+ if debug:
+ print(f"Dict: Interface{ind}.gli_term[{neighbour}]['wetting']=\n",interface.gli_term[neighbour]['wetting'])
+ print("should be the same as \n")
+ print(f"From Dict Interface{ind}.gli_term_prev[{subdomain}]['wetting']=\n",interface.gli_term_prev[subdomain]['wetting'])
gli_readout_tmp = np.zeros(gli_assembled_tmp[phase].size, dtype = float)
interface.write_gli_dofs(to_array = gli_readout_tmp, #
interface_dofs = interface_dofs[phase],#
@@ -493,6 +528,8 @@ class DomainPatch(df.SubDomain):
# this is set to True because we need gli_prev
previous_iter = True
)
+ if debug:
+ print(f"read gli_prev_iter[{subdomain}]:\n", gli_readout_tmp)
# read neighbouring pressure values from interface and write them to
# the function self.neighbouring_interface_pressure
@@ -511,7 +548,7 @@ class DomainPatch(df.SubDomain):
gli_pressure_part = df.assemble(-2*dt*Lambda[phase]*pa_prev*phi_a*ds)
gli_p_part = gli_pressure_part.get_local()
gli_prev_neighbour = gli_readout_tmp
- gli = gli_p_part + gli_prev_neighbour
+ gli = gli_p_part -gli_prev_neighbour
# check if there are shared dofs with another interface.
if dofs_in_common_with_other_interfaces[phase].size != 0:
for common_dof in dofs_in_common_with_other_interfaces[phase]:
@@ -532,7 +569,8 @@ class DomainPatch(df.SubDomain):
interface.gli_term_prev[subdomain].update(#
{'nonwetting': interface.gli_term[neighbour]['nonwetting'].copy()}
)
- else:
+
+ if neighbour_iter_num + 1 < iteration:
print(f"neighbour_iter_num = {neighbour_iter_num} for neighbour = {neighbour}\n",
f"and iteration = {iteration} on subdomain {subdomain}\n",
"It is expected that this case should only be happening if the \n",
@@ -553,6 +591,9 @@ class DomainPatch(df.SubDomain):
for ind in self.has_interface:
# self._calc_gli_term_on(interface, iteration)
interface = self.interface[ind]
+ if debug:
+ print(f"Interface{ind}.gli_term[{subdomain}][{phase}]=\n",interface.gli_term[subdomain][phase])
+ print(f"before writing to interface[{ind}] gli_assembled[{phase}]=\n",self.gli_assembled[phase])
# write gli dofs to interface dicts for communiction
interface.read_gli_dofs(from_array = self.gli_assembled[phase], #
interface_dofs = interface_dofs[phase],#
@@ -562,6 +603,8 @@ class DomainPatch(df.SubDomain):
subdomain_ind = subdomain,#
previous_iter = False
)
+ if debug:
+ print(f"gli_terms after writing to interface[{ind}] Interface[{ind}].gli_term[{subdomain}][{phase}]=\n",interface.gli_term[subdomain][phase])
### END calc_gli_term
def null_all_interface_iteration_numbers(self) -> None:
@@ -611,7 +654,7 @@ class DomainPatch(df.SubDomain):
'nonwetting' : df.interpolate(df.Constant(0), self.function_space['nonwetting'])
}
- def _init_dof_and_vertex_maps(self) -> None:
+ def _init_dof_maps(self) -> None:
""" calculate dof maps and the vertex to parent map.
This method sets the class Variables
@@ -622,15 +665,18 @@ class DomainPatch(df.SubDomain):
mesh_data = self.mesh.data()
# local to parent vertex index map
self.parent_mesh_index = mesh_data.array('parent_vertex_indices',0)
+ print(f"on subdomain{self.subdomain_index} parent_vertex_indices: {self.parent_mesh_index}")
self.dof2vertex = dict()
self.vertex2dof = dict()
for phase in self.has_phases:
self.dof2vertex.update(#
{phase :df.dof_to_vertex_map(self.function_space[phase])}#
)
+ print(f"on subdomain{self.subdomain_index} dof2vertex: {self.dof2vertex[phase]}")
self.vertex2dof.update(#
{phase :df.vertex_to_dof_map(self.function_space[phase])}#
)
+ print(f"on subdomain{self.subdomain_index} vertex2dof: {self.vertex2dof[phase]}")
def _init_markers(self) -> None:
""" define boundary markers
@@ -645,7 +691,9 @@ class DomainPatch(df.SubDomain):
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:
- marker_value = glob_index + 1
+ # 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)
# create outerBoundary objects and mark outer boundaries with 1
@@ -661,7 +709,11 @@ class DomainPatch(df.SubDomain):
internal = False,
tol= self.tol) #self.mesh.hmin()/100
)
- self.outer_boundary[index].mark(self.outer_boundary_marker, 1)
+ # similar to the interfaces self.outer_boundary_def_points is a
+ # dictionary enumerating all the boundary parts. The enumeration
+ # starts at 0, so we set the marker value to index+1
+ boundary_marker_value = index+1
+ self.outer_boundary[index].mark(self.outer_boundary_marker, boundary_marker_value)
def _init_measures(self):
@@ -682,7 +734,7 @@ class DomainPatch(df.SubDomain):
marker = self.interface_marker
for ind in self.has_interface:
# the marker value for the interfaces is always global interface index+1
- marker_value = ind + 1
+ marker_value = self.interface[ind].marker_value
self._dof_indices_of_interface.update(#
{ind: dict()}
)
@@ -730,6 +782,70 @@ class DomainPatch(df.SubDomain):
)
+ def write_solution_to_xdmf(self, file: str, #
+ time: float = None,#
+ write_iter_for_fixed_time: bool = False,#
+ ):
+ """
+ Weither write the solution of the simulation and corresponding time
+ or the solution and corresponding iteration at fixed time
+ (write_iter_for_fixed_time = True) to disk.
+ This is needed to visualize the solutions.
+
+ PARAMETERS
+ file str File name of file to write to
+ time float time at which to write the solution
+ 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:
+ self.pressure[phase].rename("pressure_"+"{}".format(phase), "pressure_"+"{}".format(phase))
+ file.write(self.pressure[phase], t)
+ else:
+ i = self.iteration_number
+ self.pressure[phase].rename("pressure_"+"{}".format(phase), #
+ "pressure(t="+"{}".format(t)+")_"+"{}".format(phase)
+ )
+ file.write(self.pressure[phase], i)
+ # create function to copy the subsequent errors to
+ error = df.Function(self.function_space[phase])
+ error.assign(self.pressure[phase] - self.pressure_prev_iter[phase])
+ error.rename("pi_"+"{}".format(phase)+"-pim1_"+"{}".format(phase), #
+ "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)
+ # # write unknowns to file
+ # rho.rename("rho","density")
+ # v.rename("v","velocity")
+ # phi.rename("phi","phasefield")
+ # mu.rename("mu","dfdc")
+ # tau.rename("tau","dfdrho")
+ # sigma.rename("sigma","gradc")
+ #
+ # for var in [rho, v, phi ,mu, tau, sigma]:
+ # file.write(var, self.t)
+ #
+ # # write pressure to file
+ # p = df.project(self.pressure(),df.FunctionSpace(self._mesh, self._element.sub_elements()[0]),solver_type="bicgstab")
+ # if self.project_to_cg == True:
+ # p = df.project(self.pressure(),df.FunctionSpace(self._mesh, self._element.sub_elements()[0].reconstruct(family="CG")),solver_type="bicgstab")
+ #
+ # p.rename("p","pressure")
+ # file.write(p, self.t)
+
+
+
# END OF CLASS DomainPatch
@@ -940,6 +1056,7 @@ class Interface(BoundaryPart):
"""
def __init__(self, #
global_index: int,#
+ # marker_value: int,#
adjacent_subdomains: np.array = None,#
# this is an array of bools
isRichards: np.array = None,#
@@ -949,6 +1066,7 @@ class Interface(BoundaryPart):
self.adjacent_subdomains = adjacent_subdomains
# global index ob the subdomain
self.global_index = global_index
+ self.marker_value = self.global_index+1
# pressure_values is a dictionary of dictionaries, one for each of the
# adjacent subdomains, holding the dof values over the interface of the
# pressure (solution) values for communication. It is initialised by the
@@ -1004,11 +1122,15 @@ class Interface(BoundaryPart):
according to self.inside
interface_marker_value: #type int marker value of interface_marker
"""
- interface_vertex_indices = self._vertex_indices(interface_marker, interface_marker_value)
+ interface_vertex_indices = self._vertex_indices(interface_marker, interface_marker_value,
+ print_vertex_indices = False)
+ print(f"local interface vertices: {interface_vertex_indices}")
dof2vertex = df.dof_to_vertex_map(function_space)
+ print(f"dof2vertex map: {dof2vertex}\n")
# check which dofs actually lie on the interface using the dof2vertex map.
is_a_dof_on_interface = np.isin(dof2vertex, interface_vertex_indices, assume_unique=True)
-
+ print(f"is_a_dof_on_interface: {is_a_dof_on_interface}")
+ print(f"this means the following dofs are interface dofs:", np.nonzero(is_a_dof_on_interface)[0])
return np.nonzero(is_a_dof_on_interface)[0]
def init_interface_values(self,#
@@ -1020,7 +1142,7 @@ class Interface(BoundaryPart):
The MeshFunction interface_marker is supposed to be one defined on the
global mesh, so that the indices used in the dictionary are the global
- vertex indices and not indicis of some submesh.
+ vertex indices and not indices of some submesh.
# Parameters
interface_marker: #type df.MeshFunction marker function marking
@@ -1030,7 +1152,7 @@ class Interface(BoundaryPart):
"""
vertex_indices = self._vertex_indices(interface_marker, #
interface_marker_value,#
- print_vertex_indices = False)
+ print_vertex_indices = True)
self.pressure_values = dict()
#self.pressure_values is a dictionay which contains for each index in
# adjacent_subdomains a dictionary for wetting and nonwetting phase with
@@ -1104,13 +1226,13 @@ class Interface(BoundaryPart):
# from_function.vector() is orderd by dof and not by vertex numbering of the mesh.
# parent needs a mesh vertex index of the submesh, therefore d2v is needed.
self.pressure_values[subdomain_ind][phase].update({parent[d2v[dof_index]]: from_function.vector()[dof_index]})
- print(f"have written subdomain{subdomain_ind} pressure to interface{self.global_index} for phase {phase}\n", self.pressure_values[subdomain_ind][phase])
+ # print(f"have written subdomain{subdomain_ind} pressure to interface{self.global_index} for phase {phase}\n", self.pressure_values[subdomain_ind][phase])
else:
for dof_index in interface_dofs:
# from_function.vector() is orderd by dof and not by vertex numbering of the mesh.
# parent needs a mesh vertex index of the submesh, therefore d2v is needed.
self.pressure_values_prev[subdomain_ind][phase].update({parent[d2v[dof_index]]: from_function.vector()[dof_index]})
- print(f"have written subdomain{subdomain_ind} pressure_prev_iter to interface{self.global_index} for phase {phase}\n", self.pressure_values_prev[subdomain_ind][phase])
+ # print(f"have written subdomain{subdomain_ind} pressure_prev_iter to interface{self.global_index} for phase {phase}\n", self.pressure_values_prev[subdomain_ind][phase])
def read_gli_dofs(self, from_array: np.array, #
@@ -1153,13 +1275,13 @@ class Interface(BoundaryPart):
# from_array is orderd by dof and not by vertex numbering of the mesh.
# parent needs a mesh vertex index of the submesh, therefore d2v is needed.
self.gli_term[subdomain_ind][phase].update({parent[d2v[dof_index]]: from_array[dof_index]})
- print(f"have written subdomain{subdomain_ind} gl dofs to interface{self.global_index} for phase {phase}\n", self.gli_term[subdomain_ind][phase])
+ # print(f"have written subdomain{subdomain_ind} gl dofs to interface{self.global_index} for phase {phase}\n", self.gli_term[subdomain_ind][phase])
else:
for dof_index in interface_dofs:
# from_array is orderd by dof and not by vertex numbering of the mesh.
# parent needs a mesh vertex index of the submesh, therefore d2v is needed.
self.gli_term_prev[subdomain_ind][phase].update({parent[d2v[dof_index]]: from_array[dof_index]})
- print(f"have written subdomain{subdomain_ind} gl_prev dofs to interface{self.global_index} for phase {phase}\n", self.gli_term_prev[subdomain_ind][phase])
+ # print(f"have written subdomain{subdomain_ind} gl_prev dofs to interface{self.global_index} for phase {phase}\n", self.gli_term_prev[subdomain_ind][phase])
def write_pressure_dofs(self, to_function: df.Function, #
@@ -1209,15 +1331,15 @@ class Interface(BoundaryPart):
# from_function.vector() is orderd by dof and not by vertex numbering of the mesh.
# parent needs a mesh vertex index of the submesh, therefore d2v is needed.
to_function.vector()[dof_index] = self.pressure_values[subdomain_ind][phase][parent[d2v[dof_index]]]
- print("read pressure interface values\n", self.pressure_values[subdomain_ind][phase])
- print("wrote these dofs to function \n", to_function.vector()[:])
+ # print("read pressure interface values\n", self.pressure_values[subdomain_ind][phase])
+ # print("wrote these dofs to function \n", to_function.vector()[:])
else:
for dof_index in interface_dofs:
# from_function.vector() is orderd by dof and not by vertex numbering of the mesh.
# parent needs a mesh vertex index of the submesh, therefore d2v is needed.
to_function.vector()[dof_index] = self.pressure_values_prev[subdomain_ind][phase][parent[d2v[dof_index]]]
- print("read previous pressure interface values\n", self.pressure_values_prev[subdomain_ind][phase])
- print("wrote these dofs to function \n", to_function.vector()[:])
+ # print("read previous pressure interface values\n", self.pressure_values_prev[subdomain_ind][phase])
+ # print("wrote these dofs to function \n", to_function.vector()[:])
def write_gli_dofs(self, to_array: np.array, #
@@ -1265,15 +1387,15 @@ class Interface(BoundaryPart):
# to_array is orderd by dof and not by vertex numbering of the mesh.
# parent needs a mesh vertex index of the submesh, therefore d2v is needed.
to_array[dof_index] = self.gli_term[subdomain_ind][phase][parent[d2v[dof_index]]]
- print("read gl interface values\n", self.gli_term[subdomain_ind][phase])
- print("wrote these dofs to array \n", to_array[:])
+ # print("read gl interface values\n", self.gli_term[subdomain_ind][phase])
+ # print("wrote these dofs to array \n", to_array[:])
else:
for dof_index in interface_dofs:
# to_array is orderd by dof and not by vertex numbering of the mesh.
# parent needs a mesh vertex index of the submesh, therefore d2v is needed.
to_array[dof_index] = self.gli_term_prev[subdomain_ind][phase][parent[d2v[dof_index]]]
- print("read gl interface values\n", self.gli_term_prev[subdomain_ind][phase])
- print("wrote these dofs to array \n", to_array[:])
+ # print("read gl interface values\n", self.gli_term_prev[subdomain_ind][phase])
+ # print("wrote these dofs to array \n", to_array[:])
#### Private Methods ####
@@ -1318,7 +1440,7 @@ class Interface(BoundaryPart):
interface_vertex_number += 1
if print_vertex_indices:
- print("mesh node indices: \n", vertex_indices)
+ print(f"mesh node indices: \n", vertex_indices)
return vertex_indices
### End Class Interface