From aed6c8e003dd67549a94debe852d9068e6447422 Mon Sep 17 00:00:00 2001
From: David Seus <david.seus@ians.uni-stuttgart.de>
Date: Wed, 3 Apr 2019 18:34:58 +0200
Subject: [PATCH] write Lsolver_step
---
LDDsimulation/LDDsimulation.py | 83 ++++++++++++++-----------
LDDsimulation/domainPatch.py | 62 +++++++++---------
RR-2-patch-test-case/RR-2-patch-test.py | 8 ++-
3 files changed, 85 insertions(+), 68 deletions(-)
diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index a0f5976..c6f9da6 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -91,18 +91,18 @@ class LDDsimulation(object):
# 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
- # turned into df.Expressions objects by the method update_DirichletBC_dictionaries()
+ # turned into df.Expressions objects by the method update_DirichletBC_dictionary()
# this is mostly in place when examples with exact solutions are calculated.
self.dirichletBC_Expressions = None
# dictionary to hold the df.Expressions for the external boundary created
# from the above self.dirichletBC_Expressions. These are
- # turned into df.dirichletBC objects by the method update_DirichletBC_dictionaries()
+ # turned into df.dirichletBC objects by the method update_DirichletBC_dictionary()
self.dirichletBC_dfExpression = dict()
# dictionary to store the df.dirichletBC objects in. These are used by the
# solver to set the boundary conditions
self.outerBC = dict()
# Flag to toggle between a case with exact solution and a purely numerical
- # simulation. This changes the behaviour of self.update_DirichletBC_dictionaries()
+ # simulation. This changes the behaviour of self.update_DirichletBC_dictionary()
self.exact_solution_case = False
# dictionary holding the source expressions for each subdomain. These are
# saved to the respective subdomains by the self._init_subdomains() method.
@@ -190,7 +190,7 @@ class LDDsimulation(object):
self._eval_sources(interpolation_degree = 2,
time = time,
subdomain_index = ind)
- self.update_DirichletBC_dictionaries()
+ 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.
@@ -287,7 +287,9 @@ class LDDsimulation(object):
# so gli_term_assembled needs to actually be subtracted from the rhs.
rhs_assembled = rhs_without_gli_assembled - gli_term_assembled
# apply outer Dirichlet boundary conditions if present.
-
+ if subdomain.outer_boundary is not None:
+ self.outerBC[subdomain_index]['wetting'].apply(form_assembled, rhs_assembled)
+ solution_pw = subdomain.pressure['wetting'].vector()
else:
print("implement two phase assembly")
@@ -504,42 +506,49 @@ class LDDsimulation(object):
"""
num = subdomain_index
subdomain = self.subdomain[num]
- mesh = subdomain.mesh
- V = subdomain.function_space
- # Here the dictionary has to be created first because t = 0.
- if np.abs(time) < self.tol :
- self.outerBC.update({num: dict()})
- self.dirichletBC_dfExpression.update({num: dict()})
-
- pDC = self.dirichletBC_Expressions[num]
- boundary_marker = subdomain.outer_boundary_marker
- if subdomain.isRichards:
- # note that the default interpolation degree is 2
+ # we only need to set Dirichlet boundary conditions if we have an outer
+ # boundary.
+ 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.tol :
- # time = 0 so the Dolfin Expression has to be created first.
- pDCw = df.Expression(pDC['wetting'], domain = mesh, degree = interpolation_degree, t=time)
- self.dirichletBC_dfExpression[num].update({'wetting': pDCw})
- else:
- pDCw = self.dirichletBC_dfExpression[num]['wetting'].t = time
+ self.outerBC.update({num: dict()})
+ self.dirichletBC_dfExpression.update({num: dict()})
- self.outerBC[num].update({'wetting': df.DirichletBC(V['wetting'], pDCw, boundary_marker, 1)})
- else:
- if np.abs(time) < self.tol :
- # time = 0 so the Dolfin Expression has to be created first.
- pDCw = df.Expression(pDC['wetting'], domain = mesh, degree = interpolation_degree, t=time)
- pDCnw = df.Expression(pDC['nonwetting'], domain = mesh, degree = interpolation_degree, t=time)
- self.dirichletBC_dfExpression[num].update(
- {'wetting': pDCw},#
- {'nonwetting': pDCnw},#
- )
+ pDC = self.dirichletBC_Expressions[num]
+ boundary_marker = subdomain.outer_boundary_marker
+ if subdomain.isRichards:
+ # 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.
+ pDCw = df.Expression(pDC['wetting'], domain = mesh, degree = interpolation_degree, t=time)
+ self.dirichletBC_dfExpression[num].update({'wetting': pDCw})
+ else:
+ pDCw = self.dirichletBC_dfExpression[num]['wetting'].t = time
+
+ self.outerBC[num].update({'wetting': df.DirichletBC(V['wetting'], pDCw, boundary_marker, 1)})
else:
- pDCw = self.dirichletBC_dfExpression[num]['wetting'].t = time
- pDCnw = self.dirichletBC_dfExpression[num]['nonwetting'].t = time
+ if np.abs(time) < self.tol :
+ # time = 0 so the Dolfin Expression has to be created first.
+ pDCw = df.Expression(pDC['wetting'], domain = mesh, degree = interpolation_degree, t=time)
+ pDCnw = df.Expression(pDC['nonwetting'], domain = mesh, degree = interpolation_degree, t=time)
+ self.dirichletBC_dfExpression[num].update(
+ {'wetting': pDCw},#
+ {'nonwetting': pDCnw},#
+ )
+ else:
+ pDCw = self.dirichletBC_dfExpression[num]['wetting'].t = time
+ pDCnw = self.dirichletBC_dfExpression[num]['nonwetting'].t = time
- self.outerBC[num].update(#
- {'wetting': df.DirichletBC(V['wetting'], pDCw, boundary_marker, 1)},#
- {'nonwetting': df.DirichletBC(V['nonwetting'], pDCnw, boundary_marker, 1)},#
- )
+ self.outerBC[num].update(#
+ {'wetting': df.DirichletBC(V['wetting'], pDCw, boundary_marker, 1)},#
+ {'nonwetting': df.DirichletBC(V['nonwetting'], pDCnw, boundary_marker, 1)},#
+ )
+ 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, interpolation_degree: int = 2, #
time: float = 0,
diff --git a/LDDsimulation/domainPatch.py b/LDDsimulation/domainPatch.py
index b4ba1d6..ac93c5e 100644
--- a/LDDsimulation/domainPatch.py
+++ b/LDDsimulation/domainPatch.py
@@ -216,12 +216,12 @@ class DomainPatch(df.SubDomain):
subdomain_ind = self.subdomain_index)
- def govering_problem(self, phase: str, iter_num: int = 0) -> tp.Dict[str, fl.Form]:
- """ return the governing form and right hand side for phase phase and iteration
- iter_num as a dictionary.
+ def govering_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 and iteration
- iter_num (without the gli terms) 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
@@ -248,17 +248,11 @@ class DomainPatch(df.SubDomain):
# we need to have all interfaces in the form
interface_forms = []
for interface in self.has_interface:
- interface_forms.append((dt*Lambda*pw*phi_w)*ds[interface])
+ interface_forms.append((dt*Lambda*pw*phi_w)*ds(interface)
form1 = (Lw*pw*phi_w)*dx
form2 = dt/mu_w*df.dot(kw(S(pw_prev_iter))*df.grad(pw), df.grad(v))*dx
form = form1 + form2 + df.sum(interface_forms)
- # # assemble rhs for
- # rhs_interface_forms = []
- # self._calc_gli_term(iteration = iter_num)
- # # gli = self.gli_assembled['wetting']
- # for index in self.has_interface:
- # gli = self.interface[index].gli_term[self.subdomain_index]['wetting']
- # rhs_interface_forms.append((dt*gli*phi_w)*ds[index])
+ # # assemble rhs
rhs1 = (Lw*pw_prev_iter*phi_w)*dx
rhs2 = -(porosity*(S(pw_prev_iter) - S(pw_prev_timestep))*phi_w)*dx
rhs_without_gli = rhs1 + rhs2 + dt*source_w*phi_w*dx
@@ -268,8 +262,8 @@ class DomainPatch(df.SubDomain):
}
return form_and_rhs
else:
- La = self.L['alpha']
- pa = self.pressure['alpha']
+ La = self.L['phase']
+ pa = self.pressure['phase']
# this is pw_{i-1} in the paper
pw_prev_iter = self.pressure_prev_iter['wetting']
pnw_prev_iter = self.pressure_prev_iter['nonwetting']
@@ -277,12 +271,12 @@ class DomainPatch(df.SubDomain):
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
- phi_a = self.testfunction['alpha']
- Lambda_a = self.lambda_param['alpha']
- ka = self.relative_permeability['alpha']
+ phi_a = self.testfunction['phase']
+ Lambda_a = self.lambda_param['phase']
+ ka = self.relative_permeability['phase']
S = self.saturation
- mu_a = self.viscosity['alpha']
- source_a = self.source['alpha']
+ mu_a = self.viscosity['phase']
+ source_a = self.source['phase']
# 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
@@ -290,7 +284,7 @@ class DomainPatch(df.SubDomain):
# gather interface forms
interface_forms = []
for interface in self.has_interface:
- interface_forms.append((dt*Lambda_a*pa*phi_a)*ds[interface])
+ interface_forms.append((dt*Lambda_a*pa*phi_a)*ds(interface))
form1 = (La*pa*phi_a)*dx
form2 = dt/mu_a*df.dot(ka(S(pc_prev_iter))*df.grad(pa), df.grad(phi_a))*dx
@@ -314,9 +308,9 @@ class DomainPatch(df.SubDomain):
# 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.append((df.Constant(0)*phi_a)*ds[interface])
+ interface_forms.append((df.Constant(0)*phi_a)*ds(interface))
else:
- interface_forms.append((dt*Lambda_a*pa*phi_a)*ds[interface])
+ interface_forms.append((dt*Lambda_a*pa*phi_a)*ds(interface))
form1 = (La*pa*phi_a)*dx
form2 = dt/mu_a*df.dot(ka(1 - S(pc_prev_iter))*df.grad(pa), df.grad(phi_a))*dx
@@ -338,10 +332,10 @@ class DomainPatch(df.SubDomain):
'Expected either phase == wetting or phase == nonwetting ')
return None
- def calc_gli_term(self, iteration: int) -> None: # tp.Dict[str, fl.Form]
- """calculate the gl terms for the iteration number iteration
+ def calc_gli_term(self) -> 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 iteration and save them
+ calculate the gl terms for the iteration number of the subdomain and save them
to self.gli_assembled. The term gets saved as an already assembled sparse
ds form to be added to the rhs by the solver. This is due to the need of
having to communicate dofs over the interface. Additionally, the dofs are
@@ -350,19 +344,20 @@ class DomainPatch(df.SubDomain):
Caution: The array self.gli_assembled needs to be added to the rhs with
a minus sign by the solver!
"""
+ iteration = self.iteration_number
+ subdomain = self.subdomain_index
+ dt = self.timestep_size
+ Lambda = self.lambda_param
+
for ind in self.has_interface:
# self._calc_gli_term_on(interface, iteration)
interface = self.interface[ind]
- subdomain = self.subdomain_index
# neighbour index
neighbour = interface.neighbour[subdomain]
neighbour_iter_num = interface.current_iteration[neighbour]
- dt = self.timestep_size
ds = self.ds(ind)
# needed for the read_gli_dofs() functions
interface_dofs = self._dof_indices_of_interface[ind],#
-
- Lambda = self.lambda_param
if iteration == 0:
n = df.FacetNormal(self.mesh)
pw_prev_iter = self.pressure_prev_timestep
@@ -585,6 +580,13 @@ class DomainPatch(df.SubDomain):
# finally we are done and can step up the current iteration count.
interface.current_iteration[subdomain] += 1
+ def null_all_interface_iteration_numbers(self) -> None:
+ """ reset interface.current_iteration[subdomain] = 0 for all interfaces
+ of the subdomain.
+ """
+ subdomain = self.subdomain_index
+ for index in self.has_interface:
+ self.interface[index].current_iteration[subdomain] = 0
#### PRIVATE METHODS
def _init_function_space(self) -> None:
diff --git a/RR-2-patch-test-case/RR-2-patch-test.py b/RR-2-patch-test-case/RR-2-patch-test.py
index 12938f2..55af1aa 100755
--- a/RR-2-patch-test-case/RR-2-patch-test.py
+++ b/RR-2-patch-test-case/RR-2-patch-test.py
@@ -67,6 +67,9 @@ subdomain_def_points = [sub_domain0_vertices,#
# in the below list, index 0 corresponds to the 12 interface which has index 1
interface_def_points = [interface12_vertices]
+# if a subdomain has no outer boundary write None instead, i.e.
+# i: None
+# if i is the index of the inner subdomain.
outer_boundary_def_points = {
# subdomain number
1 : subdomain1_outer_boundary_verts,
@@ -167,6 +170,8 @@ exact_solution = {
2: {'wetting': '1.0 - (1.0 + t*t)*(1.0 + x[1]*x[1])'}
}
+# similary to the outer boundary dictionary, if a patch has no outer boundary
+# None should be written instead of an expression.
dirichletBC = exact_solution
# def saturation(pressure, subdomain_index):
@@ -195,7 +200,8 @@ simulation.set_parameters(output_dir = "",#
timestep_size = timestep_size,#
sources = source_expression,#
initial_conditions = initial_condition,#
- outer_dirichletBC = dirichletBC,#
+ dirichletBC_Expressions = dirichletBC,#
+ exact_solution_case = True,#
)
simulation.initialise()
--
GitLab