Skip to content
Snippets Groups Projects
Commit aed6c8e0 authored by David Seus's avatar David Seus
Browse files

write Lsolver_step

parent 2af00011
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
LDDsimulation/LDDsimulation.py
+ 46
− 37
View file @ aed6c8e0
......@@ -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,6 +506,9 @@ class LDDsimulation(object):
"""
num = subdomain_index
subdomain = self.subdomain[num]
# 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.
......@@ -540,6 +545,10 @@ class LDDsimulation(object):
{'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,
......
LDDsimulation/domainPatch.py
+ 32
− 30
View file @ aed6c8e0
......@@ -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:
......
RR-2-patch-test-case/RR-2-patch-test.py
+ 7
− 1
View file @ aed6c8e0
......@@ -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()
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment