Select Git revision
domainPatch.py
domainPatch.py 64.17 KiB
"""Class definitions for subdomains representing soil patches for LDD-TPR simulation
This file contains class definitions for the class DomainPatch which defines the
subdomains used for the LDD simulation. Each subdomain of this class holds
information about permeability, porosity, the model/physics used on the patch
as well as information about which parts of the boundary are outer boundary of
the global domain and which parts are interfaces to neighboring soil patches of
the same class.
Additionally, a class BoundaryPart is defined, which accepts a List of dolfin
Points vertices, representing the vertices of a 1D polygone. The class contains
the method BoundaryPart.inside(meshpoint,on_boundary) which is used to mark
faces of a given subdomain mesh that lie on the 1D-polygone defined by vertices.
Assumed is the use of python 3.6
"""
import dolfin as df
import mshr
import numpy as np
# import domainPatch as dp
import typing as tp
import ufl as fl
import LDDsimulation
import boundary_and_interface as bi
from termcolor import colored
from solutionFile import SolutionFile
# Interface = tp.NewType('Interface', tp.any)
class DomainPatch(df.SubDomain):
""" DomainPatch Class defines subdomains used in the LDD simulation
The class DomainPatch defines subdomains for the LDD simulation. Each subdomain
of this class holds information about the phyisics involved like viscosity,
porosity, the model/physics that is assumed on the patch
as well as information about which parts of the boundary are outer boundary of
the global domain and which parts are interfaces to neighboring soil patches of
the same class.
The constructor accepts
# Parameters
isRichards #type: bool return true if Richards model is
assumed on patch.
mesh #type df::Mesh & Dolfin mesh object, submesh of the
subdomain.
porosity #type float porosity of the the soil in the patch
assumed to be constant in space.
viscosity #type tp.List[float] viscosity(ies) of the the
fluid phases in the patch viscosity[0] is assumed
to be the viscosity of the wetting phase and if
present (isRichards == False) viscosity[1] is
assumed to be the viscosity of the non-wetting phase.
interfaces: #type tp.List[bi.Interface] List of Interface class
objects from the LDDsimulation
class.
has_interface #type tp.List[int]: list of indices (w.r.t. the global
numbering of interfaces in
LDDsimulation.interfaces) indicating
which interfaces are part of the
subdomain.
L #type tp.List[float] L parameters for the L-scheme. In
case isRichards == False (i.e
two-phase flow is assumed on the
patch) we have two parameters, one
for each equation.
lambda_param #tp.List[float] L parameters for the L-scheme. In
case isRichards == False (i.e
two-phase flow is assumed on the
patch) we have two parameters, one
for each equation.
densities tp.Dict[int, tp.Dict[str, float]] Dictionary for both phases
holding the densities for the gravity
term.
include_gravity bool flag to toggle wether or not to
include the gravity term
gravity_acceleration float gravity acceleration for the gravity
term
## Public Variables
isRichards #type bool set by input, see above
mesh #type df.Mesh set by input, see above
porosity #type float set by input, see above
viscosity #type tp.List[float] set by input, see above
interfaces: #type tp.List[bi.Interface] set by input, see above
has_interface #type tp.List[int] set by input, see above
L #type tp.List[float] set by input, see above
lambda_param #type tp.List[float] set by input, see above
function_space #type tp.Dict[str: df.Function] function space
used for wetting
and nonwetting
pressures. This is
set by
self._init_function_space()
parent_mesh_index #type np.array parent vertex map of the submesh.
boundary_marker #type df.MeshFunction marker function to mark all
boundaries of the subdomain.
The marker value marking
interface i with global index
has_interface[i] will be
has_interface[i] aswell.
## Public Methods
"""
### constructor
def __init__(self, #
subdomain_index: int,#
isRichards: bool,#
mesh: df.Mesh,#
porosity: float,#
viscosity: tp.Dict[str, float],#
outer_boundary_def_points: tp.Dict[str, tp.List[df.Point]],#
# tp.List[Interface] doesn't work, because it hasen't been defined yet.
interfaces: tp.List[tp.Type[bi.Interface]],#
has_interface: tp.List[int],#
L: tp.Dict[str, float],#
lambda_param: tp.Dict[str, float],#
relative_permeability: tp.Dict[str, tp.Callable[...,None]],#
saturation: tp.Callable[..., None],#
timestep_size: float,#
densities: tp.Dict[str, float] = None,
include_gravity: bool = False,#
gravity_acceleration: float = 9.81,
interpolation_degree: int = 10,
tol: float = None,#
degree: int = 1,
):
# because we declare our own __init__ method for the derived class BoundaryPart,
# we overwrite the __init__-method of the parent class df.SubDomain. However,
# we need the methods from the parent class, so we call the __init__-method
# from df.SubDomain to access it.
df.SubDomain.__init__(self)
if tol:
self.tol = tol
else:
self.tol = df.DOLFIN_EPS
### Class Variables/Objects set by the input to the constructor
self.subdomain_index = subdomain_index
self.isRichards = isRichards
self.mesh = mesh
self.porosity = porosity
self.viscosity = viscosity
self.outer_boundary_def_points = outer_boundary_def_points
self.interface = interfaces
self.has_interface = has_interface
self.densities = densities
self.include_gravity = include_gravity
self.gravity_acceleration = gravity_acceleration
self.L = L
self.lambda_param = lambda_param
self.relative_permeability = relative_permeability
self.saturation = saturation
# timestep size, tau in the paper
self.timestep_size = timestep_size
self.interpolation_degree = interpolation_degree
self.FEM_Lagrange_degree = degree
### 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 = ['nonwetting', 'wetting'] #['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()
# dictionary holding for each interface in self.has_interface the dof
# indices that are shared with another interface (possibly none). This info is needed for
# the calculation of the gli terms.
self._interface_has_common_dof_indices = dict()
# List of objects of clas outerBoundary initialised by self._init_markers()
# can be NONE if no outer boundary is present.
self.outer_boundary = []
# measures are set by self._init_measures()
self.iteration_number = 0
self.dx = None
self.ds = None
# solution function(s) for the form set by LDDsimulation._init_initial_values()
self.pressure = dict()
# this is being populated by LDDsimulation._eval_sources()
self.source = dict()
# dictionary of FEM function for each phase holding the previous iteration.
self.pressure_prev_iter = dict()
# dictionary of FEM function for each phase holding the pressures of the previous timestep t_n.
self.pressure_prev_timestep = dict()
# test function(s) for the form set by _init_function_space
self.testfunction = None
# Dictionary holding the exact solution expression if we have a case with exact solution.
# Is set by LDDsimualtion._init_exact_solution_expression()
self.pressure_exact = None
# Dictionary of FEM function of self.function_space["pressure"] holding the interface
# dofs of the previous iteration of the neighbouring pressure. This is used
# to assemble the gli terms in the method self.calc_gli_term().
# self.neighbouring_interface_pressure = None
# valiable holding a solver object of type df.KrylovSolver. It is set the
# first time LDDsimulation.LDDsolver is run.
self.linear_solver = None
self._init_function_space()
self._init_dof_maps()
self._init_markers()
self._init_measures()
if self.include_gravity:
self.gravity_term = dict()
self._calc_gravity_expressions()
else:
self.gravity_term = None
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
# MAGIC METHODS
def __str__(self):
"""print out a few usefull informations of the subdomain"""
model = "two-phase"
if self.isRichards:
model = "Richards"
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"
else:
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."\
+ outer_boundary_str\
+ interface_string
return string
#### PUBLIC METHODS
def write_pressure_to_interfaces(self, previous_iter: int = False):
""" save the interface values of self.pressure to all neighbouring interfaces
This method is used by the LDDsimulation.LDDsolver method.
# parameters
previous_iter bool Flag to toggle wether to write self.pressure
to interface.pressure_values or to
interface.pressure_values_prev
"""
subdomain = self.subdomain_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]})
else:
pass
def governing_problem(self, phase: str) -> tp.Dict[str, fl.Form]:
""" return the governing form and right hand side for phase phase as a dictionary.
return the governing form ant right hand side for phase phase (without the gli terms) as a dictionary.
CAUTION: the gli terms are not assembled in this method and need to be
calculated by the solver separately and subtracted from the rhs!
"""
# define measures
dx = self.dx
ds = self.ds
dt = self.timestep_size
porosity = self.porosity
if self.isRichards:
# this assumes that atmospheric pressure has been normalised to 0.
# and since
pc_prev_iter = -self.pressure_prev_iter['wetting']
pc_prev_timestep = -self.pressure_prev_timestep['wetting']
if phase == 'nonwetting':
print("CAUTION: You invoked domainPatch.governing_problem() with\n", #
"Parameter phase = 'nonwetting'. But isRichards = True for",
"current subdomain. \nResetting phase = 'wetting'.\n")
phase = 'wetting'
else:
pw_prev_iter = self.pressure_prev_iter['wetting']
pnw_prev_iter = self.pressure_prev_iter['nonwetting']
pw_prev_timestep = self.pressure_prev_timestep['wetting']
pnw_prev_timestep = self.pressure_prev_timestep['nonwetting']
pc_prev_iter = pnw_prev_iter - pw_prev_iter
pc_prev_timestep = pnw_prev_timestep - pw_prev_timestep
La = self.L[phase]
pa = self.trialfunction[phase]
# this is pw_{i-1} in the paper
pa_prev_iter = self.pressure_prev_iter[phase]
phi_a = self.testfunction[phase]
Lambda_a = self.lambda_param[phase]
ka = self.relative_permeability[phase]
S = self.saturation
mu_a = self.viscosity[phase]
source_a = self.source[phase]
if self.include_gravity:
z_a = self.gravity_term[phase]
# the first part of the form and rhs are the same for both wetting and nonwetting
form1 = (La*pa*phi_a)*dx
rhs1 = (La*pa_prev_iter*phi_a)*dx
# the forms of wetting and nonwetting phase differ by a minus sign
# and by interface_form terms if the neighbour of the current patch
# assumes Richards model
if phase == "wetting":
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
if self.include_gravity:
# z_a included the minus sign already, see self._calc_gravity_expressions
rhs_gravity = dt/mu_a*df.dot(ka(S(pc_prev_iter))*df.grad(z_a), df.grad(phi_a))*dx
elif phase == "nonwetting":
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, 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
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:
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
}
return form_and_rhs
def calc_gl0_term(self, debug: bool = False) -> None: # tp.Dict[str, fl.Form]
"""calculate the gl0 terms for all interfaces of the subdomain.
calculate the gl0 terms for all interfaces of the subdomain, interpolate
them to our functionspace and save them
to the dictionries of the interfaces. This method gets called by the LDDsolver
before entering the L-iterations in each time step.
This method assumes that self.pressure_prev_timestep has assigned the values
of self.pressure for all phases.
"""
subdomain = self.subdomain_index
dt = self.timestep_size
Lambda = self.lambda_param
for interf_ind in self.has_interface:
interface = self.interface[interf_ind]
if self.isRichards:
# assuming that we normalised atmospheric pressure to zero,
# pc = -pw in the Richards case.
pc_prev = -self.pressure_prev_timestep['wetting']
else:
pw_prev = self.pressure_prev_timestep['wetting']
pnw_prev = self.pressure_prev_timestep['nonwetting']
pc_prev = pnw_prev - pw_prev
# n = df.FacetNormal(self.mesh)
S = self.saturation
for phase in self.has_phases:
# viscosity
mu_a = self.viscosity[phase]
# relative permeability
ka = self.relative_permeability[phase]
# previous timestep pressure
pa_prev = self.pressure_prev_timestep[phase]
# testfunction
# phi_a = self.testfunction[phase]
if self.include_gravity:
z_a = self.gravity_term[phase]
if phase == 'nonwetting' and interface.neighbour_isRichards[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. We don't
# need to do anything, since gli_term_prev has been
# initialised to zero.
pass
# gli_assembled_tmp[phase] += df.assemble(df.Constant(0)*ds).get_local()
else:
if phase == 'wetting':
# the wetting phase is always present and we always need
# to calculate a gl0 term.
# TODO: add intrinsic permeabilty!!!!
if self.include_gravity:
# 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)
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]
# # 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]
global_facet_ind = local2global_facet[local_facet_ind]
for gli_dof_index, gli_dof_coord in gli_dofs_along_facet.items():
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] \
- Lambda[phase]*pa_prev.vector()[p_dof_index]
if debug:
print(f"gl0.vector()[gli_dof_ind] = {gl0.vector()[gli_dof_index]}")
# save this newly calculated dof to the interface
# dictionary gli_term_prev
dof_coord_tupel = (gli_dof_coord[0], gli_dof_coord[1])
interface.gli_term_prev[subdomain][phase][global_facet_ind].update(
{dof_coord_tupel: gl0.vector()[gli_dof_index]}
)
### END calc_gl0_term
def calc_gli_term(self,
interface_index: int,
phase: str,
debug: bool = False) -> df.Function: # tp.Dict[str, fl.Form]
"""calculate the gli term for phase phase, interface interface_index
and the iteration number of the subdomain
calculate the gl term for phase phase, interface interface_index and
the iteration number of the subdomain and return it as a dolfin Function
used for assembling the governing form in self.govering_problem. Additionally,
save the newly computed gli_term to the gli dictionries of the interface.
PARAMETERS
------------------------------------
interface_index int interface index for which to calculate the
gli term.
phase str phase to calculate the gli term for.
# DEBUG: bool toggle the debug output.
IMPLEMENTATION DETAILS
Each interface has several dictionaries to store the gli terms and interface
pressures. Let interface be the interface with the index interface_index.
For each phase in self.has_phases we have four dictionaries
interface.gli_term[subdomain][phase],
interface.gli_term_prev[subdomain][phase],
interface.gli_term[neighbour][phase],
interface.gli_term_prev[neighbour][phase],
storing the gli term of the current iteration as well as the previous one,
for both the subdomain and its neighbour that it shares the interface with.
The same holds true for the pressure terms.
This is needed because the neighbouring subdomain (relative to the subdomain
we are currently on) might still need the gli term and the pressures of
the previous iteration depending on wether or not it is one iteration ahead or not.
So, depending on the current iteration number self.iteration_number (which
is updated by the solver method LDDsimulation.LDDsolver) in relation to
the iteration number of the neighbouring subdomain, we choose different
dictionaries to read the neighbouring previous pressures from.
Similarly the dictionaries used to save the newly computed gli terms must
be chosen appropriately.
The paradigm for reading the neighbouring gli_terms is, that the correct
terms we need to read are ALWAYS stored in interface.gli_term_prev.
The self.calc_gl0_method e.g. stores the gl0 values in the interface.gli_term_prev
dictionaries.
Baring in mind, that for each iteration the solver loops over the subdomains,
and the first thing it does is raise the iteration number of the subdomain
it tackles, we know that the iteration numbers can only differ by at most 1.
We get the follwing two cases:
a) subdomain_iteration == neighbour_iter_num + 1:
subdomain and neighbour started with the same iteration number when the
solver entered the new iteration. The solver raised the subodomain iteration
number by 1, so our subdomain is now ahead of the neighbouring subdomain.
The following steps are carried out.
1. step: read previous neighbouring pressure values from interface
dictionary. Since in this case the neighbour has not yet iterated beyond
the iteration we are currently in,
pressure_prev = interface.pressure_values[neighbour][phase]
holds the "previous pressure" we need.
2. step: choose the dictionary to save the gli_term in once they are
calculated. Since in this case, the neighbour has not done the next
iteration the neighbouring subdomain will need the values stored in
interface.gli_term_prev[subdomain][phase]
of the current subdomain. Therefore we cannot overwrite these values
and save the newly calculated gli_terms to
gli_save_dictionary = interface.gli_term[subdomain][phase]
b) subdomain_iteration == neighbour_iter_num:
1. step: read previous neighbouring pressure values from interface
dictionary. In this case the neighbour has iterated once beyond the
iteration we are currently in, so
pressure_prev = interface.pressure_values_prev[neighbour][phase]
holds the "previous pressure", since the LDDsolver always writes
newly calculated pressures to the dictionary interface.pressure_values.
2. step: choose the dictionary to save the gli_term in once they are
calculated.
In this case, the neighbour has done the next iteration already
will has already used the previous gli term stored in gli_term_prev
of the current subdomain. Therefor, we save the the newly calculated
gli terms to
gli_save_dictionary = interface.gli_term_prev[subdomain][phase]
After calculating gli:
Now, since our neighbour previously was in situation a) it saved the
gli_term to interface.gli_term[neighbour] and this will be overwritten
in the next step, if we don't save it to gli_term_prev of the neighbour.
Similarly for the pressure values. So in the present case, we need to
save interface.gli_term[neighbour][phase] to
interface.gli_term_prev[neighbour][phase]
and interface.pressure_values[neighbour][phase] to
interface.pressure_values_prev[neighbour][phase]
"""
# this is the number of the current iteration of the subdomain.
iteration = self.iteration_number
subdomain = self.subdomain_index
interface = self.interface[interface_index]
# gli should be initialised by zero.
gli = df.Function(self.function_space["gli"][phase])
# neighbour index
neighbour = interface.neighbour[subdomain]
# update the current_iteration number of subdomain stored in the interface
# to the current iteration number of the subdomain.
interface.current_iteration[subdomain] = iteration
# save the iteration number of the neighbour
neighbour_iter_num = interface.current_iteration[neighbour]
if (phase == 'wetting') or (not interface.neighbour_isRichards[subdomain]):
# read previous neighbouring pressure values from interface.
# depending on our iteration number and the iteration number of
# the neighbour, weed to either read from the previous pressure
# dictionary or the current one.
if iteration == neighbour_iter_num + 1:
# in this case the neighbour has not yet iterated beyond
# the iteration we are currently in, so
# interface.pressure_values[neighbour] holds the "previous pressure"
# we need to read these values from the interface to calculate the current gli.
pressure_prev = interface.pressure_values[neighbour][phase]
# choose dictionary to save the newly calculated dofs to.
# in this case, the neighbour has not done the next iteration
# and will need both pressure and gli_term_prev of the current
# subdomain. Therefor, we save the the calculated gli_terms to
# the gli_term dictionary of the current subdomain, not the
# gli_term_prev dictionary.
gli_save_dictionary = interface.gli_term[subdomain][phase]
elif iteration == neighbour_iter_num:
# this is the only possible other case. In this case the
# neighbour has iterated once beyond the iteration we are
# currently in, so interface.pressure_values_prev holds the "previous pressure"
# we need to read from the interface to calculate gli.
pressure_prev = interface.pressure_values_prev[neighbour][phase]
# choose dictionary to save the newly calculated dofs to.
# This should be the case iteration == neighbour_iter_num:
# in this case, the neighbour has done the next iteration already
# and will not need both pressure_prev and gli_term_prev of the current
# subdomain. Therefor, we save the the calculated gli_terms to
# the gli_term_prev dictionary of the current subdomain, not the
# gli_term dictionary.
gli_save_dictionary = interface.gli_term_prev[subdomain][phase]
else:
raise RuntimeError("\n!!!!!!!!! Warning !!!!!!!!!!!!!! \n"+\
"this case should not appear\n"+"neighbour_iter_num =" + "{}".format(neighbour_iter_num)\
+ "for neighbour = "+ "{}".format(neighbour)+\
"and iteration = " + "{} on subdomain".format(iteration)+ "{}".format(subdomain)\
+ "You suck! Revisite DomainPatch.calc_gli_term() and LDDsimulation.LDDsolver ASAP.")
# 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]
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]
# read gli_prev term from the neighbour.
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.
# 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]])
# 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
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 = 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]])
# 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
gli_dof_coord_tupel = (gli_dof_coord[0], gli_dof_coord[1])
gli_save_dict.update({gli_dof_coord_tupel: gli_value})
# In the following case, the gli_term of the neighbour is saved to
# gli_term currently and will be overwritten in the next step,
# if we don't save it to gli_term_prev of the neighbour.
# similarly for the pressure values.
if iteration == neighbour_iter_num:
# print(f"we are on subdomain{subdomain} and interface{interface.global_index} (alt index:{ind}) has neighbouring \n",
# f"subdomains: subdomain{subdomain} and subdomain{neighbour}")
interface.gli_term_prev[neighbour][phase].update(#
{global_facet_index: interface.gli_term[neighbour][phase][global_facet_index].copy()}
)
interface.pressure_values_prev[neighbour][phase].update(
{global_facet_index: interface.pressure_values[neighbour][phase][global_facet_index].copy()}
)
return gli
### END calc_gli_term
def null_all_interface_iteration_numbers(self) -> None:
""" reset interface.current_iteration[subdomain] = 0 for all interfaces
of the subdomain.
"""
subdomain = self.subdomain_index
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, degree: int = None) -> None:
""" create function space for solution and trial functions
This method sets the class variable
self.function_space["pressure"]
self.trialfunction
self.pressure
self.testfunction
INPUT
degree int degree of the Lagrange ansatz space
"""
if degree is None:
degree = self.FEM_Lagrange_degree
function_name = ["pressure", "gli", "flux"]
self.function_space = dict()
for function in function_name:
self.function_space.update({function: dict()})
if function == "pressure":
self.trialfunction = dict()
self.testfunction = dict()
# we want the same space for both primary variables (pw
# and pnw) to be able to calculate pnw-pw without issues.
pressure_space = df.FunctionSpace(self.mesh, 'P', degree)
for phase in self.has_phases:
if function == "pressure":
self.function_space[function].update({phase: pressure_space})
self.trialfunction.update({phase: df.TrialFunction(self.function_space["pressure"][phase])})
self.testfunction.update({phase: df.TestFunction(self.function_space["pressure"][phase])})
elif function == "gli":
degree = self.function_space["pressure"][phase].ufl_element().degree()
self.function_space[function].update({phase: df.FunctionSpace(self.mesh, 'DG', degree)})
else:
# here function = flux
degree = self.function_space["pressure"][phase].ufl_element().degree()
self.function_space[function].update({phase: df.VectorFunctionSpace(self.mesh, 'DG', degree)})
def _init_dof_maps(self) -> None:
""" calculate dof maps and the vertex to parent map.
This method sets the class Variables
self.dof2vertex
self.vertex2dof
"""
mesh_data = self.mesh.data()
# # 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()
self.dofmap = dict()
self.dofmap.update({"pressure": dict()})
self.dofmap.update({"gli": dict()})
self.dofmap.update({"flux": dict()})
for phase in self.has_phases:
self.dof2vertex.update(#
{phase :df.dof_to_vertex_map(self.function_space["pressure"][phase])}#
)
# print(f"on subdomain{self.subdomain_index} dof2vertex: {self.dof2vertex[phase]}")
self.vertex2dof.update(#
{phase :df.vertex_to_dof_map(self.function_space["pressure"][phase])}#
)
# print(f"on subdomain{self.subdomain_index} vertex2dof: {self.vertex2dof[phase]}")
self.dofmap["pressure"].update(
{phase: self.function_space["pressure"][phase].dofmap()}
)
self.dofmap["gli"].update(
{phase: self.function_space["gli"][phase].dofmap()}
)
self.dofmap["flux"].update(
{phase: dict()}
)
self.dofmap["flux"][phase].update(
{"x": self.function_space["flux"][phase].sub(0).dofmap()}
)
self.dofmap["flux"][phase].update(
{"y": self.function_space["flux"][phase].sub(1).dofmap()}
)
def _init_markers(self) -> None:
""" define boundary markers
This method sets the class Variables
self.interface_marker
self.outer_boundary_marker
"""
self.outer_boundary_marker = df.MeshFunction('size_t', self.mesh, self.mesh.topology().dim()-1)
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)
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.
if self.outer_boundary_def_points is None:
self.outer_boundary = None
else:
for index, boundary_points in self.outer_boundary_def_points.items():
self.outer_boundary.append(#
bi.BoundaryPart(vertices=boundary_points,#
internal=False,
tol=self.tol,
marker_value=index+1
) #self.mesh.hmin()/100
)
# 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
print(f"on subdomain{self.subdomain_index} I have boundary part {index}")
boundary_marker_value = self.outer_boundary[index].marker_value
self.outer_boundary[index].mark(
self.outer_boundary_marker, boundary_marker_value
)
def _init_measures(self):
""" define measures for the governing form
This method sets the class Variables
self.dx
self.ds
"""
# domain measure
self.dx = df.Measure('dx', domain = self.mesh)
# measure over the interfaces
self.ds = df.Measure('ds', domain = self.mesh, subdomain_data = self.interface_marker)
def _calc_interface_dof_indices_and_coordinates(self, debug: bool = True):
""" calculate dictionaries containing for each local facet index of
interfaces a dictionary containing {interface_dof_index: dof_coordinates}
"""
function_name = ["pressure", "gli", "flux"]
marker = self.interface_marker
self.interface_dof_indices_and_coordinates = dict()
for function in function_name:
self.interface_dof_indices_and_coordinates.update(
{function: dict()}
)
function_space = self.function_space[function]
for interface_index in self.has_interface:
interface = self.interface[interface_index]
marker_value = interface.marker_value
self.interface_dof_indices_and_coordinates[function].update(
{interface_index: dict()}
)
if debug:
print(f"\nOn subdomain {self.subdomain_index}")
for phase in self.has_phases:
if function == "flux":
self.interface_dof_indices_and_coordinates[function][interface_index].update(#
{phase: dict()}
)
# print(f"\ndofs on interface for function space {function} for phase {phase}:")
self.interface_dof_indices_and_coordinates[function][interface_index][phase].update(#
{"x": interface.dofs_and_coordinates(function_space[phase].sub(0), marker, marker_value, debug=False)}
)
self.interface_dof_indices_and_coordinates[function][interface_index][phase].update(#
{"y": interface.dofs_and_coordinates(function_space[phase].sub(1), marker, marker_value, debug=False)}
)
else:
if debug:
print(f"\ndofs on interface for function space {function} for phase {phase}:")
self.interface_dof_indices_and_coordinates[function][interface_index].update(#
{phase: interface.dofs_and_coordinates(function_space[phase], marker, marker_value, debug=debug)}
)
def _calc_corresponding_dof_indices(self, debug=True):
""" 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 dof indices and coordinates along
# the interface. These dictionaries have the indices of
# facets belonging to the interface as keys (with respect to
# the submesh numbering) and the dictionary
# containing pairs {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 subdomain submesh.
# Global facet index refers to the facet numbering
# of the mesh of the whole computational domain.
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_y": None,
"pressure": None}}
)
# flux_x dofs
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}
)
# defensive sanity check if the dof has been set.
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")
# flux_y dofs
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}
)
# sanity check for flux_y
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")
# pressure dofs
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")
if debug:
# routine to check if dofs were correctly matched. Here, we compare
# the coordinates corresponding to the different indices. If they
# match, we correctly identified the dof indices.
for interf_ind in self.has_interface:
interface = self.interface[interf_ind]
for phase in self.has_phases:
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]
for local_facet_ind in interface.outer_normals[subdomain].keys():
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 gli_dof_index, gli_dof_coord in gli_dofs_along_facet.items():
corres_flux_x_dof_ind = self.interface_corresponding_dof_index[interf_ind][phase][local_facet_ind][gli_dof_index]["flux_x"]
corres_flux_y_dof_ind = self.interface_corresponding_dof_index[interf_ind][phase][local_facet_ind][gli_dof_index]["flux_y"]
corres_pressure_dof_ind = self.interface_corresponding_dof_index[interf_ind][phase][local_facet_ind][gli_dof_index]["pressure"]
corres_flux_x_dof_coord = flux_x_dofs_along_facet[corres_flux_x_dof_ind]
corres_flux_y_dof_coord = flux_y_dofs_along_facet[corres_flux_y_dof_ind]
corres_pressure_dof_coord = p_dofs_along_facet[corres_pressure_dof_ind]
flux_x_a_matches_b = np.allclose(corres_flux_x_dof_coord ,gli_dof_coord, rtol=1e-10, atol=1e-14)
flux_x_b_matches_a = np.allclose(gli_dof_coord, corres_flux_x_dof_coord, rtol=1e-10, atol=1e-14)
if flux_x_a_matches_b and flux_x_b_matches_a:
print(f"flux_x: gli dof ({gli_dof_index}, {gli_dof_coord}) correctly matched to flux_x dof ({corres_flux_x_dof_ind}, {corres_flux_x_dof_coord})")
else:
print(f"in flux_x")
raise RuntimeError(f"np.allclose(a,b) != np.allclose(b,a) while comparing ({corres_flux_x_dof_coord}) and ({gli_dof_coord})",
"dofs matching might be faulty.")
flux_y_a_matches_b = np.allclose(corres_flux_y_dof_coord ,gli_dof_coord, rtol=1e-10, atol=1e-14)
flux_y_b_matches_a = np.allclose(gli_dof_coord, corres_flux_y_dof_coord, rtol=1e-10, atol=1e-14)
if flux_y_a_matches_b and flux_y_b_matches_a:
print(f"flux_y: gli dof ({gli_dof_index}, {gli_dof_coord}) correctly matched to flux_x dof ({corres_flux_y_dof_ind}, {corres_flux_y_dof_coord})")
else:
print(f"in flux_y")
raise RuntimeError(f"np.allclose(a,b) != np.allclose(b,a) while comparing ({corres_flux_y_dof_coord}) and ({gli_dof_coord})",
"dofs matching might be faulty.")
pressure_a_matches_b = np.allclose(corres_pressure_dof_coord ,gli_dof_coord, rtol=1e-10, atol=1e-14)
pressure_b_matches_a = np.allclose(gli_dof_coord, corres_pressure_dof_coord, rtol=1e-10, atol=1e-14)
if pressure_a_matches_b and pressure_b_matches_a:
print(f"pressure: gli dof ({gli_dof_index}, {gli_dof_coord}) correctly matched to flux_x dof ({corres_pressure_dof_ind}, {corres_pressure_dof_coord})")
else:
print(f"in pressure")
raise RuntimeError(f"np.allclose(a,b) != np.allclose(b,a) while comparing ({corres_pressure_dof_coord}) and ({gli_dof_coord})",
"dofs matching might be faulty.")
def _init_interface_values(self):
""" allocate the dictionaries for each interface used to store the
interface values for communication.
"""
# we don't need to communicate flux values directly.
marker = self.interface_marker
for interf in self.has_interface:
marker_value = self.interface[interf].marker_value
space = self.function_space
self.interface[interf].init_interface_values(
interface_marker=marker,
interface_marker_value=marker_value,
function_space=space,
subdomain_index=self.subdomain_index,
has_phases=self.has_phases,
)
def _calc_gravity_expressions(self):
""" calculate the gravity term if it is included"""
g = df.Constant(self.gravity_acceleration) #
for phase in self.has_phases:
rho = df.Constant(self.densities[phase])
self.gravity_term.update(
{phase: df.Expression('-g*rho*x[1]', domain = self.mesh, #
degree = self.interpolation_degree,
g = g, rho = rho)}
)
def write_solution_to_xdmf(self, file: tp.Type[SolutionFile], #
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)
"""
t = 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().set_local(pc_temp)
pc_temp = self.pressure["nonwetting"] - self.pressure["wetting"]
capillary_pressure.assign(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)
)
file.write(self.pressure[phase], i)
# create function to copy the subsequent errors to
error = df.Function(self.function_space["pressure"][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)
# END OF CLASS DomainPatch