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

test q1 implementation

parent 86623bae
No related branches found
No related tags found
No related merge requests found
...@@ -588,32 +588,19 @@ class LDDsimulation(object): ...@@ -588,32 +588,19 @@ class LDDsimulation(object):
""" """
subdomain = self.subdomain[subdomain_index] subdomain = self.subdomain[subdomain_index]
iteration = subdomain.iteration_number iteration = subdomain.iteration_number
# calculate the gli terms on the right hand side and store
subdomain.calc_gli_term()
for phase in subdomain.has_phases: for phase in subdomain.has_phases:
# extract L-scheme form and rhs (without gli term) from subdomain. # extract L-scheme form and rhs (without gli term) from subdomain.
governing_problem = subdomain.governing_problem(phase = phase) governing_problem = subdomain.governing_problem(phase = phase)
form = governing_problem['form'] form = governing_problem['form']
rhs_without_gli = governing_problem['rhs_without_gli'] rhs = governing_problem['rhs']
# assemble the form and rhs # assemble the form and rhs
form_assembled = df.assemble(form) form_assembled = df.assemble(form)
rhs_without_gli_assembled = df.assemble(rhs_without_gli) rhs_assembled = df.assemble(rhs)
# access the assembled gli term for the rhs.
gli_term_assembled = subdomain.gli_assembled[phase]
# subdomain.calc_gli_term() asslembles gli but on the left hand side
# so gli_term_assembled needs to actually be subtracted from the rhs.
if debug and subdomain.mesh.num_cells() < 36:
print("\nSystem before applying outer boundary conditions:")
print(f"phase = {phase}: rhs_without_gli:\n", rhs_without_gli_assembled.get_local())
# we need the rhs as a dolfin function, so we overwrite
# rhs_without_gli_assembled with the gli_term subtracted from itself.
rhs_without_gli_assembled.set_local(rhs_without_gli_assembled.get_local() - gli_term_assembled)
rhs_assembled = rhs_without_gli_assembled
if debug and subdomain.mesh.num_cells() < 36: if debug and subdomain.mesh.num_cells() < 36:
# print(f"phase = {phase}: form_assembled:\n", form_assembled.array()) # print(f"phase = {phase}: form_assembled:\n", form_assembled.array())
# print(f"phase = {phase}: rhs_without_gli:\n", rhs_without_gli_assembled.get_local()) # print(f"phase = {phase}: rhs:\n", rhs_assembled.get_local())
print(f"phase = {phase}: gli_term:\n", gli_term_assembled) print(f"phase = {phase}: gli_term:\n", gli_term_assembled)
print(f"phase = {phase}: rhs_assembled:\n", rhs_assembled.get_local()) print(f"phase = {phase}: rhs_assembled:\n", rhs_assembled.get_local())
...@@ -627,12 +614,6 @@ class LDDsimulation(object): ...@@ -627,12 +614,6 @@ class LDDsimulation(object):
print("\nSystem after applying outer boundary conditions:") print("\nSystem after applying outer boundary conditions:")
# print(f"phase = {phase}: form_assembled:\n", form_assembled.array()) # print(f"phase = {phase}: form_assembled:\n", form_assembled.array())
print(f"phase = {phase}: rhs_assembled:\n", rhs_assembled.get_local()) print(f"phase = {phase}: rhs_assembled:\n", rhs_assembled.get_local())
# # set previous iteration as initial guess for the linear solver
# pi_prev = subdomain.pressure_prev_iter[phase].vector().get_local()
# subdomain.pressure[phase].vector().set_local(pi_prev)
if debug and subdomain.mesh.num_cells() < 36:
print("\npressure before solver:\n", subdomain.pressure[phase].vector().get_local()) print("\npressure before solver:\n", subdomain.pressure[phase].vector().get_local())
linear_solver = subdomain.linear_solver linear_solver = subdomain.linear_solver
...@@ -644,6 +625,7 @@ class LDDsimulation(object): ...@@ -644,6 +625,7 @@ class LDDsimulation(object):
if debug and subdomain.mesh.num_cells() < 36: if debug and subdomain.mesh.num_cells() < 36:
print("\npressure after solver:\n", subdomain.pressure[phase].vector().get_local()) print("\npressure after solver:\n", subdomain.pressure[phase].vector().get_local())
def _init_solution_files(self): def _init_solution_files(self):
""" set up solution files for saving the solution of the LDD scheme for """ set up solution files for saving the solution of the LDD scheme for
each subdomain. each subdomain.
......
This diff is collapsed.
...@@ -11,41 +11,42 @@ import functools as ft ...@@ -11,41 +11,42 @@ import functools as ft
##### Domain and Interface #### ##### Domain and Interface ####
# global simulation domain domain # global simulation domain domain
sub_domain0_vertices = [df.Point(0.0,-1.0), # sub_domain0_vertices = [df.Point(-1.0,-1.0), #
df.Point(1.0,-1.0),# df.Point(1.0,-1.0),#
df.Point(1.0,1.0),# df.Point(1.0,1.0),#
df.Point(0.0,1.0)] df.Point(-1.0,1.0)]
# interface between subdomain1 and subdomain2 # interface between subdomain1 and subdomain2
interface12_vertices = [df.Point(0.0, 0.0), interface12_vertices = [df.Point(-1.0, 0.0),
df.Point(1.0, 0.0) ] df.Point(1.0, 0.0) ]
# subdomain1. # subdomain1.
sub_domain1_vertices = [interface12_vertices[0], sub_domain1_vertices = [interface12_vertices[0],
interface12_vertices[1], interface12_vertices[1],
df.Point(1.0,1.0), sub_domain0_vertices[2],
df.Point(0.0,1.0) ] sub_domain0_vertices[3] ]
# vertex coordinates of the outer boundaries. If it can not be specified as a # vertex coordinates of the outer boundaries. If it can not be specified as a
# polygon, use an entry per boundary polygon. This information is used for defining # polygon, use an entry per boundary polygon. This information is used for defining
# the Dirichlet boundary conditions. If a domain is completely internal, the # the Dirichlet boundary conditions. If a domain is completely internal, the
# dictionary entry should be 0: None # dictionary entry should be 0: None
subdomain1_outer_boundary_verts = { subdomain1_outer_boundary_verts = {
0: [interface12_vertices[1], # 0: [interface12_vertices[1],
df.Point(1.0,1.0), # sub_domain0_vertices[2],
df.Point(0.0,1.0), # sub_domain0_vertices[3], #
interface12_vertices[0]] interface12_vertices[0]]
} }
# subdomain2 # subdomain2
sub_domain2_vertices = [df.Point(0.0,-1.0), sub_domain2_vertices = [sub_domain0_vertices[0],
df.Point(1.0,-1.0), sub_domain0_vertices[1],
interface12_vertices[1], interface12_vertices[1],
interface12_vertices[0] ] interface12_vertices[0] ]
subdomain2_outer_boundary_verts = { subdomain2_outer_boundary_verts = {
0: [interface12_vertices[0], # 0: [interface12_vertices[0], #
df.Point(0.0,-1.0), # sub_domain0_vertices[0],
df.Point(1.0,-1.0), # sub_domain0_vertices[1],
interface12_vertices[1]] interface12_vertices[1]]
} }
# subdomain2_outer_boundary_verts = { # subdomain2_outer_boundary_verts = {
# 0: [interface12_vertices[0], df.Point(0.0,0.0)],# # 0: [interface12_vertices[0], df.Point(0.0,0.0)],#
# 1: [df.Point(0.0,0.0), df.Point(1.0,0.0)], # # 1: [df.Point(0.0,0.0), df.Point(1.0,0.0)], #
...@@ -85,10 +86,10 @@ isRichards = { ...@@ -85,10 +86,10 @@ isRichards = {
} }
##################### MESH ##################################################### ##################### MESH #####################################################
mesh_resolution = 5 mesh_resolution = 10
##################### TIME ##################################################### ##################### TIME #####################################################
timestep_size = 0.01 timestep_size = 0.01
number_of_timesteps = 300 number_of_timesteps = 100
# decide how many timesteps you want analysed. Analysed means, that we write out # decide how many timesteps you want analysed. Analysed means, that we write out
# subsequent errors of the L-iteration within the timestep. # subsequent errors of the L-iteration within the timestep.
number_of_timesteps_to_analyse = 11 number_of_timesteps_to_analyse = 11
...@@ -107,6 +108,7 @@ densities = {# ...@@ -107,6 +108,7 @@ densities = {#
2 : {'wetting' :1} 2 : {'wetting' :1}
} }
# gravity_acceleration = 9.81
porosity = {# porosity = {#
# subdom_num : porosity # subdom_num : porosity
...@@ -122,8 +124,8 @@ L = {# ...@@ -122,8 +124,8 @@ L = {#
lambda_param = {# lambda_param = {#
# subdom_num : lambda parameter for the L-scheme # subdom_num : lambda parameter for the L-scheme
1 : {'wetting' :40},# 1 : {'wetting' :4},#
2 : {'wetting' :40} 2 : {'wetting' :4}
} }
## relative permeabilty functions on subdomain 1 ## relative permeabilty functions on subdomain 1
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment