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

finish writing domainPatch.governing_problem

parent 8a2e7011
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Showing
with 184 additions and 115 deletions
domainPatch.py
+ 184
− 115
View file @ 0e89dd3e
...@@ -196,13 +196,17 @@ class DomainPatch(df.SubDomain): ...@@ -196,13 +196,17 @@ class DomainPatch(df.SubDomain):
ds = self.ds ds = self.ds
dt = self.timestep_size dt = self.timestep_size
if self.isRichards: if self.isRichards:
if phase == 'nonwetting':
print("CAUTION: You invoked domainPatch.governing_problem() with \n#
Parameter phase = 'nonwetting'. But isRichards = True for \n#
current subdomain. \nReturning wetting phase equation only.\n")
Lw = self.L['wetting'] Lw = self.L['wetting']
# this is pw_i in the paper # this is pw_i in the paper
pw = self.pressure['wetting'] pw = self.pressure['wetting']
# this is pw_{i-1} in the paper # this is pw_{i-1} in the paper
pw_prev_iter = self.pressure_prev_iter['wetting'] pw_prev_iter = self.pressure_prev_iter['wetting']
pw_prev_timestep = self.pressure_prev_timestep['wetting'] pw_prev_timestep = self.pressure_prev_timestep['wetting']
v = self.testfunction['wetting'] phi_w = self.testfunction['wetting']
Lambda = self.lambda_param['wetting'] Lambda = self.lambda_param['wetting']
kw = self.relative_permeability['wetting'] kw = self.relative_permeability['wetting']
S = self.saturation S = self.saturation
...@@ -210,30 +214,101 @@ class DomainPatch(df.SubDomain): ...@@ -210,30 +214,101 @@ class DomainPatch(df.SubDomain):
source_w = self.source['wetting'] source_w = self.source['wetting']
# we need to have all interfaces in the form # we need to have all interfaces in the form
interface_forms = [] interface_forms = []
for index in self.has_interface: for interface in self.has_interface:
interface_forms.append((dt*Lambda*pw*v)*ds[index]) interface_forms.append((dt*Lambda*pw*phi_w)*ds[interface])
form = (Lw*pw*v)*dx + (dt/mu_w*df.dot(kw(S(pw_prev_iter))*df.grad(pw), df.grad(v)))*dx + df.sum(interface_forms) form1 = (Lw*pw*phi_w)*dx
# assemble rhs for form2 = dt/mu_w*df.dot(kw(S(pw_prev_iter))*df.grad(pw), df.grad(v))*dx
rhs_interface_forms = [] form = form1 + form2 + df.sum(interface_forms)
self._calc_gli_term(iteration = iter_num) # # assemble rhs for
# gli = self.gli['wetting'] # rhs_interface_forms = []
for index in self.has_interface: # self._calc_gli_term(iteration = iter_num)
gli = self.interface[index].gli_term[self.subdomain_index]['wetting'] # # gli = self.gli['wetting']
rhs_interface_forms.append((dt*gli*v)*ds[index]) # for index in self.has_interface:
rhs = (Lw*pw_prev_iter*v)*dx - ((S(pw_prev_iter) - S(pw_prev_timestep)*v)*dx + dt*source_w*v)*dx - df.sum(rhs_interface_forms) # gli = self.interface[index].gli_term[self.subdomain_index]['wetting']
# rhs_interface_forms.append((dt*gli*phi_w)*ds[index])
rhs1 = (Lw*pw_prev_iter*phi_w)*dx
rhs2 = -((S(pw_prev_iter) - S(pw_prev_timestep))*phi_w)*dx
rhs_without_gli = rhs1 + rhs2 + dt*source_w*phi_w*dx
form_and_rhs = {# form_and_rhs = {#
'form': form,# 'form': form,#
'rhs' : rhs 'rhs_without_gli' : rhs_without_gli
} }
return form_and_rhs
else: 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']
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
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']
if phase == "wetting": if phase == "wetting":
print("governing form wetting phase") # gather interface forms
interface_forms = []
for interface in self.has_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
form = form1 + form2 + df.sum(interface_forms)
rhs1 = (La*pw_prev_iter*phi_a)*dx
rhs2 = -((S(pc_prev_iter) - S(pc_prev_timestep))*phi_a)*dx
rhs_without_gli = rhs1 + rhs2 + dt*source_a*phi_a)*dx
form_and_rhs = {#
'form': form,#
'rhs_without_gli' : rhs_without_gli
}
return form_and_rhs
elif phase == "nonwetting": elif phase == "nonwetting":
print("governing form nonwetting phase") # gather interface forms
interface_forms = []
for interface in self.has_interface:
if self.interface[interface].neighbour_is_Richards(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.append((df.Constant(0)*phi_a)*ds[interface])
else: else:
raise RuntimeError('missmatch for input parameter phase. \nExpected either phase == wetting or phase == nonwetting ') 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
form = form1 + form2 + df.sum(interface_forms)
rhs1 = (La*pnw_prev_iter*phi_a)*dx
# the non wetting phase has a + before instead of a - before the bracket
rhs2 = ((S(pc_prev_iter) - S(pc_prev_timestep))*phi_a)*dx
rhs_without_gli = rhs1 + rhs2 + dt*source_a*phi_a)*dx
form_and_rhs = {#
'form': form,#
'rhs_without_gli' : rhs_without_gli
}
return form_and_rhs
return form else:
raise RuntimeError('missmatch for input parameter phase. \n#
Expected either phase == wetting or phase == nonwetting ')
#### PRIVATE METHODS #### PRIVATE METHODS
def _init_function_space(self) -> None: def _init_function_space(self) -> None:
...@@ -357,67 +432,67 @@ class DomainPatch(df.SubDomain): ...@@ -357,67 +432,67 @@ class DomainPatch(df.SubDomain):
'nonwetting': self.interface[ind].dofs_on_interface(V['nonwetting'], marker, ind)} 'nonwetting': self.interface[ind].dofs_on_interface(V['nonwetting'], marker, ind)}
) )
# def _calc_gli_term(self, iteration: int) -> None: # tp.Dict[str, fl.Form] 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
#
# calculate the gl terms for the iteration number iteration and save them calculate the gl terms for the iteration number iteration and save them
# to self.gli. to self.gli.
# """ """
# for ind in self.has_interface: for ind in self.has_interface:
# # self._calc_gli_term_on(interface, iteration) # self._calc_gli_term_on(interface, iteration)
# interface = self.interface[ind] interface = self.interface[ind]
# subdomain = self.subdomain_index subdomain = self.subdomain_index
# # neighbour index # neighbour index
# neighbour = interface.neighbour[subdomain] neighbour = interface.neighbour[subdomain]
# neighbour_iter_num = interface.current_iteration[neighbour] neighbour_iter_num = interface.current_iteration[neighbour]
#
# Lambda = self.lambda_param Lambda = self.lambda_param
# if iteration == 0: if iteration == 0:
# n = df.FacetNormal(self.mesh) n = df.FacetNormal(self.mesh)
# pw_prev_iter = self.pressure_prev_timestep pw_prev_iter = self.pressure_prev_timestep
# k = self.relative_permeability k = self.relative_permeability
# S = self.saturation S = self.saturation
# mu = self.viscosity mu = self.viscosity
#
# if not neighbour_iter_num == 0: if not neighbour_iter_num == 0:
# # save the gl term from the neighbour first to use it in the # save the gl term from the neighbour first to use it in the
# # next iteration # next iteration
# if self.isRichards: if self.isRichards:
# interface.gli_term_prev[subdomain].update(# interface.gli_term_prev[subdomain].update(#
# {'wetting': interface.gli_term[neighbour]['wetting']} {'wetting': interface.gli_term[neighbour]['wetting']}
# ) )
# else: else:
# interface.gli_term_prev[subdomain].update(# interface.gli_term_prev[subdomain].update(#
# {'wetting': interface.gli_term[neighbour]['wetting']},# {'wetting': interface.gli_term[neighbour]['wetting']},#
# {'nonwetting': interface.gli_term[neighbour]['nonwetting']} {'nonwetting': interface.gli_term[neighbour]['nonwetting']}
# ) )
#
# if self.isRichards: if self.isRichards:
# mu_w = mu['wetting'] mu_w = mu['wetting']
# kw = k['wetting'] kw = k['wetting']
# pwn_prev = pw_prev_iter['wetting'] pwn_prev = pw_prev_iter['wetting']
# #calc gl0 from the flux term #calc gl0 from the flux term
# flux = -1/mu_w*df.dot(kw(S(pwn_prev))*df.grad(pwn_prev) flux = -1/mu_w*kw(S(pwn_prev))*df.grad(pwn_prev)
# gl0 = (df.dot(flux, n) - Lambda['wetting']*pwn_prev) gl0 = (df.dot(flux, n) - Lambda['wetting']*pwn_prev)
# interface.gli_term[subdomain].update({'wetting': gl0}) interface.gli_term[subdomain].update({'wetting': gl0})
# else: else:
# mu_w = mu['wetting'] mu_w = mu['wetting']
# mu_nw = mu['nonwetting'] mu_nw = mu['nonwetting']
# pwn_prev = pw_prev_iter['wetting'] pwn_prev = pw_prev_iter['wetting']
# pnwn_prev = pw_prev_iter['nonwetting'] pnwn_prev = pw_prev_iter['nonwetting']
# kw = k['wetting'] kw = k['wetting']
# knw = k['nonwetting'] knw = k['nonwetting']
# #calc gl0 from the flux term #calc gl0 from the flux term
# fluxw = -1/mu_w*df.dot(kw(S(pnwn_prev - pwn_prev))*df.grad(pwn_prev) fluxw = -1/mu_w*kw(S(pnwn_prev - pwn_prev))*df.grad(pwn_prev)
# fluxnw = -1/mu_nw*df.dot(knw(1-S(pnwn_prev - pwn_prev))*df.grad(pnwn_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) gwl0 = (df.dot(fluxw, n) - Lambda['wetting']*pwn_prev)
# gnwl0 = (df.dot(fluxnw, n) - Lambda['nonwetting']*pnwn_prev) gnwl0 = (df.dot(fluxnw, n) - Lambda['nonwetting']*pnwn_prev)
# interface.gli_term[subdomain].update({'wetting': gwl0}) interface.gli_term[subdomain].update({'wetting': gwl0})
# interface.gli_term[subdomain].update({'nonwetting': gnwl0}) interface.gli_term[subdomain].update({'nonwetting': gnwl0})
#
# interface.current_iteration[subdomain] += 1 interface.current_iteration[subdomain] += 1
#
# else: # else: # iteration > 0
# if neighbour_iter_num == iteration: # if neighbour_iter_num == iteration:
# if self.isRichards: # if self.isRichards:
# interface.gli_term_prev[subdomain].update(# # interface.gli_term_prev[subdomain].update(#
...@@ -426,7 +501,7 @@ class DomainPatch(df.SubDomain): ...@@ -426,7 +501,7 @@ class DomainPatch(df.SubDomain):
# #
# pwn_prev = pw_prev_iter['wetting'] # pwn_prev = pw_prev_iter['wetting']
# #calc gl0 from the flux term # #calc gl0 from the flux term
# flux = -1/mu_w*df.dot(kw(S(pwn_prev))*df.grad(pwn_prev) # flux = -1/mu_w*kw(S(pwn_prev))*df.grad(pwn_prev)
# gl0 = (df.dot(flux, n) - Lambda['wetting']*pwn_prev) # gl0 = (df.dot(flux, n) - Lambda['wetting']*pwn_prev)
# interface.gli_term[subdomain].update({'wetting': gl0}) # interface.gli_term[subdomain].update({'wetting': gl0})
# else: # else:
...@@ -437,23 +512,17 @@ class DomainPatch(df.SubDomain): ...@@ -437,23 +512,17 @@ class DomainPatch(df.SubDomain):
# kw = k['wetting'] # kw = k['wetting']
# knw = k['nonwetting'] # knw = k['nonwetting']
# #calc gl0 from the flux term # #calc gl0 from the flux term
# fluxw = -1/mu_w*df.dot(kw(S(pnwn_prev - pwn_prev))*df.grad(pwn_prev) # fluxw = -1/mu_w*kw(S(pnwn_prev - pwn_prev))*df.grad(pwn_prev)
# fluxnw = -1/mu_nw*df.dot(knw(1-S(pnwn_prev - pwn_prev))*df.grad(pnwn_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) # gwl0 = (df.dot(fluxw, n) - Lambda['wetting']*pwn_prev)
# gnwl0 = (df.dot(fluxnw, n) - Lambda['nonwetting']*pnwn_prev) # gnwl0 = (df.dot(fluxnw, n) - Lambda['nonwetting']*pnwn_prev)
# interface.gli_term[subdomain].update({'wetting': gwl0}) # interface.gli_term[subdomain].update({'wetting': gwl0})
# interface.gli_term[subdomain].update({'nonwetting': gnwl0}) # interface.gli_term[subdomain].update({'nonwetting': gnwl0})
# #
# else: # else:
# print('') # print('uiae')
#
# interface.current_iteration[subdomain] += 1 # interface.current_iteration[subdomain] += 1
# # def _calc_gli_term_on(self, interface: int, iteration: int) -> None:
# # """calculate the gl terms for the iteration number iteration
# #
# # calculate the gl terms for the iteration number iteration and save them
# # to self.gli.
# # """
# # END is_Richards
# END OF CLASS DomainPatch # END OF CLASS DomainPatch
... ...
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please to comment