From 6d1059a3721376001d647c02b2ab2068a910b8a2 Mon Sep 17 00:00:00 2001
From: David Seus <david.seus@ians.uni-stuttgart.de>
Date: Fri, 29 Mar 2019 17:22:15 +0100
Subject: [PATCH] rewrite calc_gli_terms. pre rewrite of read_dofs and
write_dofs
---
LDDsimulation.py | 5 +-
RR-2-patch-test.py | 89 +++++++++-------
domainPatch.py | 255 +++++++++++++++++++++++++--------------------
3 files changed, 199 insertions(+), 150 deletions(-)
diff --git a/LDDsimulation.py b/LDDsimulation.py
index bc3064f..f137716 100644
--- a/LDDsimulation.py
+++ b/LDDsimulation.py
@@ -93,7 +93,7 @@ class LDDsimulation(object):
output_dir: str,#
subdomain_def_points: tp.List[tp.List[df.Point]],#
# this is actually an array of bools
- isRichards: bool,#
+ isRichards: tp.Dict[int, bool],#
interface_def_points: tp.List[tp.List[df.Point]],#
outer_boundary_def_points: tp.Dict[int, tp.Dict[int, tp.List[df.Point]]],#
adjacent_subdomains: tp.List[np.ndarray],#
@@ -267,11 +267,10 @@ class LDDsimulation(object):
def _init_subdomains(self) -> None:
""" generate list of opjects of class DomainPatch.
"""
- is_richards = self.isRichards
# The list is_richards has the same length as the nuber of subdomains.
# Therefor it is used in the subdomain initialisation loop.
- for subdom_num, isR in enumerate(is_richards, 1):
+ for subdom_num, isR in self.isRichards.items():
interface_list = self._get_interfaces(subdom_num)
self.subdomain.update(#
{subdom_num : dp.DomainPatch(#
diff --git a/RR-2-patch-test.py b/RR-2-patch-test.py
index 8fbb509..165394a 100755
--- a/RR-2-patch-test.py
+++ b/RR-2-patch-test.py
@@ -76,7 +76,10 @@ outer_boundary_def_points = {
# adjacent_subdomains[i] contains the indices of the subdomains sharing the
# interface i (i.e. given by interface_def_points[i]).
adjacent_subdomains = [[1,2]]
-isRichards = [True, True]
+isRichards = {
+ 1: True, #
+ 2: True
+ }
Tmax = 2
dt1 = 0.1
@@ -172,6 +175,7 @@ dirichletBC = exact_solution
#
# sa
+mesh_resolution = 3
# initialise LDD simulation class
simulation = ldd.LDDsimulation()
@@ -181,7 +185,7 @@ simulation.set_parameters(output_dir = "",#
interface_def_points = interface_def_points,#
outer_boundary_def_points = outer_boundary_def_points,#
adjacent_subdomains = adjacent_subdomains,#
- mesh_resolution = 2,#
+ mesh_resolution = mesh_resolution,#
viscosity = viscosity,#
porosity = porosity,#
L = L,#
@@ -251,7 +255,7 @@ coordinates1 = mesh_subdomain[1].coordinates()
coordinates2 = mesh_subdomain[2].coordinates()
submesh1_data = mesh_subdomain[1].data()
dom1_parent_vertex_indices = submesh1_data.array('parent_vertex_indices',0)
-print("map of dom1 vertices to parent vertices on global domain:\n", dom1_parent_vertex_indices)
+# print("map of dom1 vertices to parent vertices on global domain:\n", dom1_parent_vertex_indices)
#
# coordinates2 = mesh_subdomain[2].coordinates()
# print("\nvertices of subdomain2: \n", coordinates2)
@@ -263,7 +267,7 @@ dom2_parent_vertex_indices = submesh2_data.array('parent_vertex_indices',0)
# print("on domain: \n")
# interface_coordinates = interface[0].coordinates(interface_marker, 1, print_coordinates = True)
# interface_coordinates = interface[0]._vertex_indices(interface_marker, 1, print_vertex_indices = True)
-print("\non subdomain2: \n")
+# print("\non subdomain2: \n")
# dom2_interface_coordinates = interface[0].coordinates(interface_marker2, 1, print_coordinates = True)
dom2_interface_def_points = interface[0]._vertex_indices(interface_marker2, 1, print_vertex_indices = True)
@@ -273,7 +277,7 @@ print("Dofs on Interface dom1", dofs_on_interface1)
V2 = df.FunctionSpace(mesh_subdomain[2], 'P', 1)
dofs_on_interface2 = interface[0].dofs_on_interface(V2, interface_marker2, 1)
-print("Dofs on Interface dom1", dofs_on_interface2)
+# print("Dofs on Interface dom1", dofs_on_interface2)
# markermesh = interface_marker2.mesh()
# functionSpaceMesh = V2.mesh()
@@ -282,12 +286,12 @@ print("Dofs on Interface dom1", dofs_on_interface2)
dom1_d2v = df.dof_to_vertex_map(V1)
dom1_v2d = df.vertex_to_dof_map(V1)
-print("dof to vertex map for V1:\n", dom1_d2v)
-print("vertex to dof map for V1:\n", dom1_v2d)
+# print("dof to vertex map for V1:\n", dom1_d2v)
+# print("vertex to dof map for V1:\n", dom1_v2d)
dom2_d2v = df.dof_to_vertex_map(V2)
dom2_v2d = df.vertex_to_dof_map(V2)
-print("dof to vertex map for V2:\n", dom2_d2v)
-print("vertex to dof map for V2:\n", dom2_v2d)
+# print("dof to vertex map for V2:\n", dom2_d2v)
+# print("vertex to dof map for V2:\n", dom2_v2d)
testf1 = df.Constant(41)
f1test = df.interpolate(testf1, V1)
@@ -303,7 +307,7 @@ f1nodal = f1test.vector().get_local()
# for i in range(len(f1nodal)):
# f1nodal[i] = i
-print("nodal values of f1nodal\n", f1nodal)
+# print("nodal values of f1nodal\n", f1nodal)
#
# print('dofs on interface \n', f1vec[dofs_on_interface])
@@ -311,24 +315,24 @@ print("nodal values of f1nodal\n", f1nodal)
# print("f1test.vector()[:]\n", f1test.vector()[:])
# print("f1test.vector()[dofs_on_interface]\n", f1test.vector()[dofs_on_interface])
-print("constant function on dom1\n", f1vec)
-
-for i in dofs_on_interface1:
- print("local dof ind %i, parent = %i", i,dom1_parent_vertex_indices[dom1_d2v[i]])
- # f1test.vector() is enumerated by dof not by vertex
- f1test.vector()[i] = i
-# this enumerates according to the dof numbering
-print("testfunction f1test with set interface values\n", f1test.vector()[:])
-print("testfunction f1test with get_local\n", f1test.vector().get_local())
+# print("constant function on dom1\n", f1vec)
+#
+# for i in dofs_on_interface1:
+# print("local dof ind %i, parent = %i", i,dom1_parent_vertex_indices[dom1_d2v[i]])
+# # f1test.vector() is enumerated by dof not by vertex
+# f1test.vector()[i] = i
+# # this enumerates according to the dof numbering
+# print("testfunction f1test with set interface values\n", f1test.vector()[:])
+# print("testfunction f1test with get_local\n", f1test.vector().get_local())
# f1test.vector()[dom1_v2d[dom1_d2v[dofs_on_interface]]] = f1vec[dofs_on_interface]
# .vector()[:] and .vector().get_local() order values according to the vertex numbering
-print("testfunction f1test with compute_vertex_values\n", f1test.compute_vertex_values())
-print("testfunction f1test with compute_vertex_values mit d2vmap\n", f1test.compute_vertex_values()[dom1_d2v[:]])
-i=0
-for x in coordinates1:
- print(f"vertex {i}: f1test[{dom1_v2d[i]}] = {f1test.vector()[dom1_v2d[i]]}\tu({x}) = {f1test(x)}")
- i += 1
-print("\non subdomain1: \n")
+# print("testfunction f1test with compute_vertex_values\n", f1test.compute_vertex_values())
+# print("testfunction f1test with compute_vertex_values mit d2vmap\n", f1test.compute_vertex_values()[dom1_d2v[:]])
+# i=0
+# for x in coordinates1:
+# print(f"vertex {i}: f1test[{dom1_v2d[i]}] = {f1test.vector()[dom1_v2d[i]]}\tu({x}) = {f1test(x)}")
+# i += 1
+# print("\non subdomain1: \n")
print("vertices of subdomain1: \n", coordinates1)
# dom1_interface_coordinates = interface[0].coordinates(interface_marker1, 1, print_coordinates = True)
dom1_interface_def_points = interface[0]._vertex_indices(interface_marker1, 1, print_vertex_indices = True)
@@ -394,8 +398,8 @@ p1_0 = df.interpolate(p1_initial, V1)
p2_0 = df.interpolate(p2_initial, V2)
# Initial boundary value for the pressure on interface
-p_gamma1 = df.DirichletBC(V1, df.Constant(0.5), interface_marker1, 1)
-p_gamma2 = df.DirichletBC(V2, df.Constant(0.1), interface_marker2, 1)
+p_gamma1 = df.DirichletBC(V1, p1_exact, interface_marker1, 1)
+p_gamma2 = df.DirichletBC(V2, p2_exact, interface_marker2, 1)
outerBC1 = df.DirichletBC(V1, p1_exact, boundary_marker1, 1)
outerBC2 = df.DirichletBC(V2, p2_exact, boundary_marker2, 1)
@@ -419,26 +423,39 @@ lambda2 = lambda1
p1_i = p1_0
p2_i = p2_0
# initial g_i
+n1 = df.FacetNormal(mesh_subdomain[1])
+flux = -df.grad(p1_0)
+gl0 = (df.dot(flux, n1) - lambda1*p1_0)*v1*ds1(1)
+print('gl0 term assembliert', df.assemble(gl0).get_local())
g1_i = -lambda1*p1_i
g2_i = -lambda1*p2_i
for iterations in range(1):
- # + relative_permeability(saturation(p1_i, 1), 1)* +
- a1 = (df.inner(df.grad(u1), df.grad(v1)))*dx_1 + (timestep*lambda1*u1*v1)*ds1(1) #timestep*
- rhs1 = (L1*p1_i*v1)*dx1 - ((saturation(p1_i, 1) - saturation(p1_0, 1))*v1)*dx1 + (timestep*(source1 - g1_i)*v1)*ds1
+ # + +
+ a1 = (L1*u1*v1)*dx1 + (df.inner(relative_permeability(saturation(p1_i, 1), 1)*df.grad(u1), df.grad(v1)))*dx_1 + (timestep*lambda1*u1*v1)*ds1(1) #timestep*
+ rhs1 = (L1*p1_i*v1)*dx1 - ((saturation(p1_i, 1) - saturation(p1_0, 1))*v1)*dx1 + (timestep*(source1 - g1_i)*v1)*ds1(1)
+ rhssplit1 = (L1*p1_i*v1)*dx1 - ((saturation(p1_i, 1) - saturation(p1_0, 1))*v1)*dx1
+ extratrem = (timestep*(source1 - g1_i)*v1)*ds1(1)
+ splitrhs = df.assemble(rhssplit1)
+ print("splitrhs: \n", splitrhs.get_local())
+ extrat_assembled = df.assemble(extratrem)
+ print("extratrem: \n", extrat_assembled.get_local())
+ added = splitrhs + extrat_assembled
+ print("both added together:\n", added.get_local())
+
A1 = df.assemble(a1)
dsform = (timestep*lambda1*u1*v1)*ds1(1)
dsform_vec = df.assemble(dsform)
b1 = df.assemble(rhs1)
- p_gamma1.apply(A1,b1)
- outerBC1.apply(A1,b1)
+ #p_gamma1.apply(A1,b1)
+ # outerBC1.apply(A1,b1)
u = df.Function(V1)
U1 = u.vector()
df.solve(A1,U1,b1)
- print("Form:\n", A1.array())
+ # print("Form:\n", A1.array())
print("RHS:\n", b1.get_local())
- print("ds term:\n", dsform_vec.array())
- print('solution:\n', U1.get_local())
+ #print("ds term:\n", dsform_vec.array())
+ # print('solution:\n', U1.get_local())
p1_i.vector()[:] = U1
interface[0].write_dofs(from_function = p1_i, #
interface_dofs = dofs_on_interface1,#
diff --git a/domainPatch.py b/domainPatch.py
index 44b7964..de36a41 100644
--- a/domainPatch.py
+++ b/domainPatch.py
@@ -150,15 +150,19 @@ class DomainPatch(df.SubDomain):
# solution function(s) for the form set by _init_function_space
self.pressure = None
self.source = None
- # dictionary holding the previous iteration.
+ # FEM function holding the previous iteration.
self.pressure_prev_iter = None
- # dictionary holding the pressures of the previous timestep t_n.
+ # FEM function holding the pressures of the previous timestep t_n.
self.pressure_prev_timestep = None
# test function(s) for the form set by _init_function_space
- self.TestFunction = None
- # # dictionary which is initiated by _init_function_space holding the
- # # the gli term needed for the form assemply.
- # self.gli = None
+ 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 = None
self._init_function_space()
self._init_dof_and_vertex_maps()
self._init_markers()
@@ -194,7 +198,7 @@ class DomainPatch(df.SubDomain):
iter_num 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.
+ iter_num (without the gli terms) as a dictionary.
"""
# define measures
dx = self.dx
@@ -227,7 +231,7 @@ class DomainPatch(df.SubDomain):
# # assemble rhs for
# rhs_interface_forms = []
# self._calc_gli_term(iteration = iter_num)
- # # gli = self.gli['wetting']
+ # # 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])
@@ -241,10 +245,7 @@ class DomainPatch(df.SubDomain):
return form_and_rhs
else:
La = self.L['alpha']
- # Lnw = self.L['nonwetting']
- # this is pw_i/pnw_i in the paper
pa = self.pressure['alpha']
- # pnw = self.pressure['nonwetting']
# this is pw_{i-1} in the paper
pw_prev_iter = self.pressure_prev_iter['wetting']
pnw_prev_iter = self.pressure_prev_iter['nonwetting']
@@ -253,17 +254,14 @@ class DomainPatch(df.SubDomain):
pc_prev_iter = pnw_prev_iter - pw_prev_iter
pc_prev_timestep = pnw_prev_timestep - pw_prev_timestep
phi_a = self.testfunction['alpha']
- # phi_nw = self.testfunction['nonwetting']
Lambda_a = self.lambda_param['alpha']
- # Lambda_nw = self.lambda_param['nonwetting']
ka = self.relative_permeability['alpha']
- # knw = self.relative_permeability['nonwetting']
S = self.saturation
mu_a = self.viscosity['alpha']
- # mu_nw = self.viscosity['nonwetting']
source_a = self.source['alpha']
- # source_nw = self.source['nonwetting']
-
+ # 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":
# gather interface forms
interface_forms = []
@@ -316,6 +314,132 @@ 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
+
+ calculate the gl terms for the iteration number iteration 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
+ being saved to an interface dictionary exactly as the pressures.
+ """
+ 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)
+
+ Lambda = self.lambda_param
+ if iteration == 0:
+ n = df.FacetNormal(self.mesh)
+ pw_prev_iter = self.pressure_prev_timestep
+ k = self.relative_permeability
+ S = self.saturation
+ mu = self.viscosity
+
+ if not neighbour_iter_num == 0:
+ # save the gl term from the neighbour first to use it in the
+ # next iteration
+ if self.isRichards:
+ interface.gli_term_prev[subdomain].update(#
+ {'wetting': interface.gli_term[neighbour]['wetting']}
+ )
+ else:
+ interface.gli_term_prev[subdomain].update(#
+ {'wetting': interface.gli_term[neighbour]['wetting']},#
+ {'nonwetting': interface.gli_term[neighbour]['nonwetting']}
+ )
+
+ if self.isRichards:
+ mu_w = mu['wetting']
+ kw = k['wetting']
+ pw_prev = pw_prev_iter['wetting']
+ phi_w = self.testfunction['wetting']
+ #calc gl0 from the flux term
+ flux = -1/mu_w*kw(S(pw_prev))*df.grad(pw_prev)
+ gl0 = (df.dot(flux, n) - Lambda['wetting']*pw_prev)
+ self.gli_assembled = {
+ 'wetting': df.assemble(dt*gl0*phi_w*ds)
+ }
+ ### TODO write dofs to interface dict
+ #interface.gli_term[subdomain].update({'wetting': gl0})
+ else:
+ phi_w = self.testfunction['wetting']
+ mu_w = mu['wetting']
+ pw_prev = pw_prev_iter['wetting']
+ pnw_prev = pw_prev_iter['nonwetting']
+ pc_prev = pnw_prev - pw_prev
+ kw = k['wetting']
+ #calc gl0 from the flux term
+ if 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
+ fluxw = -1/mu_w*kw(S(pc_prev))*df.grad(pw_prev)
+ gwl0 = (df.dot(fluxw, n) - Lambda['wetting']*pw_prev)
+ self.gli_assembled = {
+ 'wetting': df.assemble(dt*gwl0*phi_w*ds),
+ # assemble the nonwetting phase gl0 term as zero functional
+ 'nonwetting': df.assemble(df.Constant(0)*ds)
+ }
+ ###TODO write gli_terms to interface
+ else:
+ # here we need to assemble a gli term for the nonwetting phase
+ phi_nw = self.testfunction['nonwetting']
+ mu_nw = mu['nonwetting']
+ knw = k['nonwetting']
+
+ fluxw = -1/mu_w*kw(S(pc_prev))*df.grad(pw_prev)
+ fluxnw = -1/mu_nw*knw(1-S(pc_prev))*df.grad(pnw_prev)
+ gwl0 = (df.dot(fluxw, n) - Lambda['wetting']*pw_prev)
+ gnwl0 = (df.dot(fluxnw, n) - Lambda['nonwetting']*pnw_prev)
+ self.gli_assembled = {
+ 'wetting': df.assemble(dt*gwl0*phi_w*ds),
+ 'nonwetting': df.assemble(dt*gnwl0*phi_nw*ds)
+ }
+ ###TODO write gli_terms to interface
+ # interface.gli_term[subdomain].update({'wetting': gwl0})
+ # interface.gli_term[subdomain].update({'nonwetting': gnwl0})
+
+ interface.current_iteration[subdomain] += 1
+
+ # else: # iteration > 0
+ # if neighbour_iter_num == iteration:
+ # if self.isRichards:
+ # interface.gli_term_prev[subdomain].update(#
+ # {'wetting': interface.gli_term[neighbour]['wetting']}
+ # )
+ #
+ # pw_prev = pw_prev_iter['wetting']
+ # #calc gl0 from the flux term
+ # flux = -1/mu_w*kw(S(pw_prev))*df.grad(pw_prev)
+ # gl0 = (df.dot(flux, n) - Lambda['wetting']*pw_prev)
+ # interface.gli_term[subdomain].update({'wetting': gl0})
+ # else:
+ # mu_w = mu['wetting']
+ # mu_nw = mu['nonwetting']
+ # pw_prev = pw_prev_iter['wetting']
+ # pnw_prev = pw_prev_iter['nonwetting']
+ # kw = k['wetting']
+ # knw = k['nonwetting']
+ # #calc gl0 from the flux term
+ # fluxw = -1/mu_w*kw(S(pnw_prev - pw_prev))*df.grad(pw_prev)
+ # fluxnw = -1/mu_nw*knw(1-S(pnw_prev - pw_prev))*df.grad(pnw_prev)
+ # gwl0 = (df.dot(fluxw, n) - Lambda['wetting']*pw_prev)
+ # gnwl0 = (df.dot(fluxnw, n) - Lambda['nonwetting']*pnw_prev)
+ # interface.gli_term[subdomain].update({'wetting': gwl0})
+ # interface.gli_term[subdomain].update({'nonwetting': gnwl0})
+ #
+ # else:
+ # print('uiae')
+ #
+ # interface.current_iteration[subdomain] += 1
+
+
#### PRIVATE METHODS
def _init_function_space(self) -> None:
""" create function space for solution and trial functions
@@ -324,13 +448,13 @@ class DomainPatch(df.SubDomain):
aswell. This method sets the class variable
self.function_space
self.pressure
- self.TestFunction
+ self.testfunction
"""
if self.isRichards:
self.function_space = {'wetting': df.FunctionSpace(self.mesh, 'P', 1)}
self.pressure = df.TrialFunction(self.function_space['wetting'])
self.testfunction = df.TestFunction(self.function_space['wetting'])
- # self.gli = {'wetting' : df.interpolate(df.Constant(0), self.function_space['wetting'])}
+ # self.gli_assembled = {'wetting' : df.interpolate(df.Constant(0), self.function_space['wetting'])}
else:
self.function_space = {#
'wetting' : df.FunctionSpace(self.mesh, 'P', 1),#
@@ -344,7 +468,7 @@ class DomainPatch(df.SubDomain):
'wetting' : df.TestFunction(self.function_space['wetting']),
'nonwetting' : df.TestFunction(self.function_space['nonwetting'])
}
- # self.gli = {#
+ # self.gli_assembled = {#
# 'wetting' : df.interpolate(df.Constant(0), self.function_space['wetting']),#
# 'nonwetting' : df.interpolate(df.Constant(0), self.function_space['nonwetting'])
# }
@@ -438,97 +562,6 @@ class DomainPatch(df.SubDomain):
'nonwetting': self.interface[ind].dofs_on_interface(V['nonwetting'], marker, ind)}
)
- def _calc_gli_term(self, iteration: int) -> None: # tp.Dict[str, fl.Form]
- """calculate the gl terms for the iteration number iteration
-
- calculate the gl terms for the iteration number iteration and save them
- to self.gli.
- """
- 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]
-
- Lambda = self.lambda_param
- if iteration == 0:
- n = df.FacetNormal(self.mesh)
- pw_prev_iter = self.pressure_prev_timestep
- k = self.relative_permeability
- S = self.saturation
- mu = self.viscosity
-
- if not neighbour_iter_num == 0:
- # save the gl term from the neighbour first to use it in the
- # next iteration
- if self.isRichards:
- interface.gli_term_prev[subdomain].update(#
- {'wetting': interface.gli_term[neighbour]['wetting']}
- )
- else:
- interface.gli_term_prev[subdomain].update(#
- {'wetting': interface.gli_term[neighbour]['wetting']},#
- {'nonwetting': interface.gli_term[neighbour]['nonwetting']}
- )
-
- if self.isRichards:
- mu_w = mu['wetting']
- kw = k['wetting']
- pwn_prev = pw_prev_iter['wetting']
- #calc gl0 from the flux term
- flux = -1/mu_w*kw(S(pwn_prev))*df.grad(pwn_prev)
- gl0 = (df.dot(flux, n) - Lambda['wetting']*pwn_prev)
- interface.gli_term[subdomain].update({'wetting': gl0})
- else:
- mu_w = mu['wetting']
- mu_nw = mu['nonwetting']
- pwn_prev = pw_prev_iter['wetting']
- pnwn_prev = pw_prev_iter['nonwetting']
- kw = k['wetting']
- knw = k['nonwetting']
- #calc gl0 from the flux term
- fluxw = -1/mu_w*kw(S(pnwn_prev - pwn_prev))*df.grad(pwn_prev)
- fluxnw = -1/mu_nw*knw(1-S(pnwn_prev - pwn_prev))*df.grad(pnwn_prev)
- gwl0 = (df.dot(fluxw, n) - Lambda['wetting']*pwn_prev)
- gnwl0 = (df.dot(fluxnw, n) - Lambda['nonwetting']*pnwn_prev)
- interface.gli_term[subdomain].update({'wetting': gwl0})
- interface.gli_term[subdomain].update({'nonwetting': gnwl0})
-
- interface.current_iteration[subdomain] += 1
-
- # else: # iteration > 0
- # if neighbour_iter_num == iteration:
- # if self.isRichards:
- # interface.gli_term_prev[subdomain].update(#
- # {'wetting': interface.gli_term[neighbour]['wetting']}
- # )
- #
- # pwn_prev = pw_prev_iter['wetting']
- # #calc gl0 from the flux term
- # flux = -1/mu_w*kw(S(pwn_prev))*df.grad(pwn_prev)
- # gl0 = (df.dot(flux, n) - Lambda['wetting']*pwn_prev)
- # interface.gli_term[subdomain].update({'wetting': gl0})
- # else:
- # mu_w = mu['wetting']
- # mu_nw = mu['nonwetting']
- # pwn_prev = pw_prev_iter['wetting']
- # pnwn_prev = pw_prev_iter['nonwetting']
- # kw = k['wetting']
- # knw = k['nonwetting']
- # #calc gl0 from the flux term
- # fluxw = -1/mu_w*kw(S(pnwn_prev - pwn_prev))*df.grad(pwn_prev)
- # fluxnw = -1/mu_nw*knw(1-S(pnwn_prev - pwn_prev))*df.grad(pnwn_prev)
- # gwl0 = (df.dot(fluxw, n) - Lambda['wetting']*pwn_prev)
- # gnwl0 = (df.dot(fluxnw, n) - Lambda['nonwetting']*pnwn_prev)
- # interface.gli_term[subdomain].update({'wetting': gwl0})
- # interface.gli_term[subdomain].update({'nonwetting': gnwl0})
- #
- # else:
- # print('uiae')
- #
- # interface.current_iteration[subdomain] += 1
# END OF CLASS DomainPatch
--
GitLab