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

before attempting to fix calculation of gli terms

parent b1af8421
Branches
Tags
No related merge requests found
...@@ -239,25 +239,25 @@ class LDDsimulation(object): ...@@ -239,25 +239,25 @@ class LDDsimulation(object):
raise RuntimeError('Parameters note set!\nYou forgott to run self.set_parameters(**kwds)') raise RuntimeError('Parameters note set!\nYou forgott to run self.set_parameters(**kwds)')
self._add_parameters_to_output_dirname() self._add_parameters_to_output_dirname()
self._init_analyse_timesteps() self._init_analyse_timesteps()
df.info(colored("initialise meshes and marker functions ...\n", "orange")) df.info(colored("initialise meshes and marker functions ...\n", "yellow"))
self._init_meshes_and_markers() self._init_meshes_and_markers()
df.info(colored("initialise interfaces ...\n", "orange")) df.info(colored("initialise interfaces ...\n", "yellow"))
self._init_interfaces() self._init_interfaces()
df.info(colored("initialise subdomains ...\n", "orange")) df.info(colored("initialise subdomains ...\n", "yellow"))
self._init_subdomains() self._init_subdomains()
df.info(colored("creating xdmf files for each subdomain ...\n", "orange")) df.info(colored("creating xdmf files for each subdomain ...\n", "yellow"))
self._init_solution_files() self._init_solution_files()
# initial values get initialised here to be able to call the solver at different # initial values get initialised here to be able to call the solver at different
# timesteps without having to call self.run(). # timesteps without having to call self.run().
df.info(colored("create initial values ...", "orange")) df.info(colored("create initial values ...", "yellow"))
self._init_initial_values() self._init_initial_values()
df.info(colored("create source Expressions ...", "orange")) df.info(colored("create source Expressions ...", "yellow"))
self._init_sources() self._init_sources()
# initialise exact solution expression dictionary if we have an exact # initialise exact solution expression dictionary if we have an exact
# solution case. # solution case.
self._init_exact_solution_expression() self._init_exact_solution_expression()
if self.exact_solution: if self.exact_solution:
df.info(colored("writing exact solution for all time steps to xdmf files ...", "orange")) df.info(colored("writing exact solution for all time steps to xdmf files ...", "yellow"))
self.write_exact_solution_to_xdmf() self.write_exact_solution_to_xdmf()
self._initialised = True self._initialised = True
...@@ -372,7 +372,6 @@ class LDDsimulation(object): ...@@ -372,7 +372,6 @@ class LDDsimulation(object):
# iterations in, used to analyse the behaviour of the iterates within # iterations in, used to analyse the behaviour of the iterates within
# the current time step. # the current time step.
solution_over_iteration_within_timestep = dict() solution_over_iteration_within_timestep = dict()
subsequent_errors_over_iteration = dict()
# before the iteration gets started all iteration numbers must be nulled # before the iteration gets started all iteration numbers must be nulled
# and the self._calc_isdone_for_subdomain dictionary must be set to False # and the self._calc_isdone_for_subdomain dictionary must be set to False
for ind, subdomain in self.subdomain.items(): for ind, subdomain in self.subdomain.items():
...@@ -385,20 +384,11 @@ class LDDsimulation(object): ...@@ -385,20 +384,11 @@ class LDDsimulation(object):
solution_over_iteration_within_timestep.update( # solution_over_iteration_within_timestep.update( #
{ind: SolutionFile(self._mpi_communicator, filename)} {ind: SolutionFile(self._mpi_communicator, filename)}
) )
subsequent_errors_over_iteration.update( {ind: dict()} )
# # we need a dictionries for each phase.
# for phase in subdomain.has_phases:
# subsequent_errors_over_iteration[ind].update( {phase: dict()} )
#
self._calc_isdone_for_subdomain.update({ind: False}) self._calc_isdone_for_subdomain.update({ind: False})
# reset all interface[has_interface].current_iteration[ind] = 0, for # reset all interface[has_interface].current_iteration[ind] = 0, for
# index has_interface in subdomain.has_interface # index has_interface in subdomain.has_interface
subdomain.null_all_interface_iteration_numbers() subdomain.null_all_interface_iteration_numbers()
# update the time of the source terms before starting the iteration.
# The sources and sinks are the same throughout the iteration.
# self._eval_sources(time = time,
# subdomain_index = ind)
self.update_DirichletBC_dictionary(time = time, # self.update_DirichletBC_dictionary(time = time, #
subdomain_index = ind, subdomain_index = ind,
) )
...@@ -409,9 +399,14 @@ class LDDsimulation(object): ...@@ -409,9 +399,14 @@ class LDDsimulation(object):
# evaluate the sources to the current time. # evaluate the sources to the current time.
subdomain.source[phase].t = time subdomain.source[phase].t = time
# set the solution of the previous timestep as new prev_timestep pressure # set the solution of the previous timestep as new prev_timestep pressure
# print(colored(f"id(subdomain{subdomain.subdomain_index}.pressure_prev_timestep[{phase}])", "red"), " =\n", f"{id(subdomain.pressure_prev_timestep[phase])}")
# print(f"subdomain{subdomain.subdomain_index}.pressure_prev_timestep[{phase}].vector()[:]", " =\n", f"{subdomain.pressure_prev_timestep[phase].vector()[:]}")
# print(f"subdomain{subdomain.subdomain_index}.pressure[{phase}].vector()[:]", " =\n", f"{subdomain.pressure[phase].vector()[:]}")
subdomain.pressure_prev_timestep[phase].assign( subdomain.pressure_prev_timestep[phase].assign(
subdomain.pressure[phase] subdomain.pressure[phase]
) )
# print(colored(f"id(subdomain{subdomain.subdomain_index}.pressure[{phase}])", "red"), " =\n", f"{id(subdomain.pressure[phase])}")
# print(f"subdomain{subdomain.subdomain_index}.pressure_prev_timestep[{phase}].vector()[:]", " =\n", f"{subdomain.pressure_prev_timestep[phase].vector()[:]}")
# All isdone flags need to be zero before we start iterating # All isdone flags need to be zero before we start iterating
self._calc_isdone_for_phase[ind].update({phase: False}) self._calc_isdone_for_phase[ind].update({phase: False})
...@@ -437,6 +432,7 @@ class LDDsimulation(object): ...@@ -437,6 +432,7 @@ class LDDsimulation(object):
self.Lsolver_step(subdomain_index = sd_index,# self.Lsolver_step(subdomain_index = sd_index,#
debug = debug, debug = debug,
) )
# bypass error checking for the first iteration.
if iteration > 0: if iteration > 0:
subsequent_iter_error = dict() subsequent_iter_error = dict()
for phase in subdomain.has_phases: for phase in subdomain.has_phases:
...@@ -446,17 +442,17 @@ class LDDsimulation(object): ...@@ -446,17 +442,17 @@ class LDDsimulation(object):
subsequent_iter_error.update( subsequent_iter_error.update(
{phase: error_calculated} {phase: error_calculated}
) )
# set flag to stop solving for that phase if error_criterion is met. # set flag to stop solving for that phase if error_criterion is met.
if subsequent_iter_error[phase] < error_criterion: if subsequent_iter_error[phase] < error_criterion:
self._calc_isdone_for_phase[sd_index].update({phase: True}) self._calc_isdone_for_phase[sd_index].update({phase: True})
if debug: if debug:
print(f"the subsequent error in iteration {iteration} for phase {phase} is {subsequent_iter_error[phase]}") print(f"t = {self.t}: subsequent error in iteration {iteration} for phase {phase} is {subsequent_iter_error[phase]}")
if analyse_timestep: if analyse_timestep:
# save the calculated L2 errors to csv file for plotting # save the calculated L2 errors to csv file for plotting
# with pgfplots. # with pgfplots
print(f"t = {self.t} the subsequent error in iteration {iteration} is {subsequent_iter_error[phase]}")
subsequent_error_filename = self.output_dir\ subsequent_error_filename = self.output_dir\
+self.output_filename_parameter_part[sd_index]\ +self.output_filename_parameter_part[sd_index]\
+"subsequent_iteration_errors" +"_at_time"+\ +"subsequent_iteration_errors" +"_at_time"+\
...@@ -467,28 +463,21 @@ class LDDsimulation(object): ...@@ -467,28 +463,21 @@ class LDDsimulation(object):
errors = subsequent_iter_error, errors = subsequent_iter_error,
time = time) time = time)
# this saves the actual difference function between subsequent
# iterations to analyse then in paraview later.
subsequent_errors_over_iteration[sd_index].update(
{ phase: error })
# both phases should be done before we mark the subdomain as # both phases should be done before we mark the subdomain as
# finished. # finished.
print("max(subsequent_iter_error.values()): ", max(subsequent_iter_error.values()))
total_subsequent_iter_error = max(subsequent_iter_error.values()) total_subsequent_iter_error = max(subsequent_iter_error.values())
# check if error criterion has been met # check if error criterion has been met
# print(f" total_subsequent_error in iter {iteration} = {total_subsequent_iter_error }\n")
if debug: if debug:
print(f" total_subsequent_error in iter {iteration} = {total_subsequent_iter_error }\n") # print(f" total_subsequent_error in iter {iteration} = {total_subsequent_iter_error }\n")
print(f"the maximum distance between any to points in a cell mesh.hmax() on subdomain{sd_index} is h={subdomain.mesh.hmax()}") print(f"the maximum distance between any to points in a cell mesh.hmax() on subdomain{sd_index} is h={subdomain.mesh.hmax()}")
print(f"the minimal distance between any to points in a cell mesh.hmin() on subdomain{sd_index} is h={subdomain.mesh.hmin()}") #print(f"the minimal distance between any to points in a cell mesh.hmin() on subdomain{sd_index} is h={subdomain.mesh.hmin()}")
print(f"the error criterion is tol = {error_criterion}") print(f"the error criterion is tol = {error_criterion}")
# wetting_done = self._calc_isdone_for_phase[sd_index]['wetting'] # wetting_done = self._calc_isdone_for_phase[sd_index]['wetting']
# nonwetting_done = self._calc_isdone_for_phase[sd_index]['nonwetting'] # nonwetting_done = self._calc_isdone_for_phase[sd_index]['nonwetting']
if total_subsequent_iter_error < error_criterion: if total_subsequent_iter_error < error_criterion:
if debug: # if debug:
print(f"subdomain{sd_index} reachd error tol = {error_criterion} in iteration {iteration}.") print(f"subdomain{sd_index} reached error = {total_subsequent_iter_error} < tol = {error_criterion} after {iteration} iterations.")
self._calc_isdone_for_subdomain[sd_index] = True self._calc_isdone_for_subdomain[sd_index] = True
# prepare next iteration # prepare next iteration
......
...@@ -21,6 +21,7 @@ import numpy as np ...@@ -21,6 +21,7 @@ import numpy as np
import typing as tp import typing as tp
import ufl as fl import ufl as fl
import LDDsimulation import LDDsimulation
from termcolor import colored
# Interface = tp.NewType('Interface', tp.any) # Interface = tp.NewType('Interface', tp.any)
# SolutionFile = tp.NewType('LDDsimulation.SolutionFile', tp.object) # SolutionFile = tp.NewType('LDDsimulation.SolutionFile', tp.object)
#import domainPatch as domainPatch #import domainPatch as domainPatch
...@@ -184,7 +185,7 @@ class DomainPatch(df.SubDomain): ...@@ -184,7 +185,7 @@ class DomainPatch(df.SubDomain):
# Dictionary of FEM function of self.function_space holding the interface # Dictionary of FEM function of self.function_space holding the interface
# dofs of the previous iteration of the neighbouring pressure. This is used # dofs of the previous iteration of the neighbouring pressure. This is used
# to assemble the gli terms in the method self.calc_gli_term(). # to assemble the gli terms in the method self.calc_gli_term().
self.neighbouring_interface_pressure = None # self.neighbouring_interface_pressure = None
# valiable holding a solver object of type df.KrylovSolver. It is set the # valiable holding a solver object of type df.KrylovSolver. It is set the
# first time LDDsimulation.LDDsolver is run. # first time LDDsimulation.LDDsolver is run.
self.linear_solver = None self.linear_solver = None
...@@ -256,8 +257,9 @@ class DomainPatch(df.SubDomain): ...@@ -256,8 +257,9 @@ class DomainPatch(df.SubDomain):
porosity = self.porosity porosity = self.porosity
if self.isRichards: if self.isRichards:
# this assumes that atmospheric pressure has been normalised to 0. # this assumes that atmospheric pressure has been normalised to 0.
pc_prev_iter = self.pressure_prev_iter['wetting'] # and since
pc_prev_timestep = self.pressure_prev_timestep['wetting'] pc_prev_iter = -self.pressure_prev_iter['wetting']
pc_prev_timestep = -self.pressure_prev_timestep['wetting']
if phase == 'nonwetting': if phase == 'nonwetting':
print("CAUTION: You invoked domainPatch.governing_problem() with\n", # print("CAUTION: You invoked domainPatch.governing_problem() with\n", #
"Parameter phase = 'nonwetting'. But isRichards = True for", "Parameter phase = 'nonwetting'. But isRichards = True for",
...@@ -365,7 +367,7 @@ class DomainPatch(df.SubDomain): ...@@ -365,7 +367,7 @@ class DomainPatch(df.SubDomain):
ds = self.ds(interface.marker_value) ds = self.ds(interface.marker_value)
# needed for the read_gli_dofs() functions # needed for the read_gli_dofs() functions
interface_dofs = self._dof_indices_of_interface[ind] interface_dofs = self._dof_indices_of_interface[ind]
# for each interface, this is a list of dofs indices belonging to another # for each interface, the following is a list of dofs indices belonging to another
# interface. The dofs corresponding to these indices must be multiplied # interface. The dofs corresponding to these indices must be multiplied
# by 1/2 to to get an average of the gli terms for dofs belonging to # by 1/2 to to get an average of the gli terms for dofs belonging to
# two different interfaces. # two different interfaces.
...@@ -406,8 +408,8 @@ class DomainPatch(df.SubDomain): ...@@ -406,8 +408,8 @@ class DomainPatch(df.SubDomain):
if self.isRichards: if self.isRichards:
# assuming that we normalised atmospheric pressure to zero, # assuming that we normalised atmospheric pressure to zero,
# pc = pw in the Richards case. # pc = -pw in the Richards case.
pc_prev = p_prev_iter['wetting'] pc_prev = -p_prev_iter['wetting']
else: else:
pw_prev = p_prev_iter['wetting'] pw_prev = p_prev_iter['wetting']
pnw_prev = p_prev_iter['nonwetting'] pnw_prev = p_prev_iter['nonwetting']
...@@ -458,8 +460,8 @@ class DomainPatch(df.SubDomain): ...@@ -458,8 +460,8 @@ class DomainPatch(df.SubDomain):
# if the neighbour of our subdomain (with index self.subdomain_index) # if the neighbour of our subdomain (with index self.subdomain_index)
# assumes a Richards model, we don't have gli terms for the # assumes a Richards model, we don't have gli terms for the
# nonwetting phase and assemble the terms as zero form. # nonwetting phase and assemble the terms as zero form.
# pass pass
gli_assembled_tmp[phase] += df.assemble(df.Constant(0)*ds).get_local() # gli_assembled_tmp[phase] += df.assemble(df.Constant(0)*ds).get_local()
else: else:
# in this case the neighbour assumes two-phase flow # in this case the neighbour assumes two-phase flow
# and we need to calculte a gl0 torm for the nonwetting # and we need to calculte a gl0 torm for the nonwetting
...@@ -498,10 +500,12 @@ class DomainPatch(df.SubDomain): ...@@ -498,10 +500,12 @@ class DomainPatch(df.SubDomain):
if debug: if debug:
print(f"in iteration {iteration} and neighbour iteration = {neighbour_iter_num}:\n") print(f"in iteration {iteration} and neighbour iteration = {neighbour_iter_num}:\n")
print(f"Interface{ind}.gli_term_prev[{subdomain}]['wetting']=\n",interface.gli_term_prev[subdomain]['wetting']) print(f"Interface{ind}.gli_term_prev[{subdomain}]['wetting']=\n",interface.gli_term_prev[subdomain]['wetting'])
interface.gli_term_prev[subdomain].update(# interface.gli_term_prev[subdomain].update(#
{'wetting': interface.gli_term[neighbour]['wetting'].copy()} {'wetting': interface.gli_term[neighbour]['wetting'].copy()}
) )
# print(colored(f"id(interface.gli_term_prev[{subdomain}][wetting])= {id(interface.gli_term_prev[subdomain]['wetting'])}","red"))
# print(colored(f"id(interface.gli_term[{neighbour}][wetting])= {id(interface.gli_term[neighbour]['wetting'])}","red"))
if debug: if debug:
print(f"Interface{ind}.gli_term[{neighbour}]['wetting']=\n",interface.gli_term[neighbour]['wetting']) print(f"Interface{ind}.gli_term[{neighbour}]['wetting']=\n",interface.gli_term[neighbour]['wetting'])
print(f"Interface{ind}.gli_term_prev[{subdomain}]['wetting']=\n",interface.gli_term_prev[subdomain]['wetting']) print(f"Interface{ind}.gli_term_prev[{subdomain}]['wetting']=\n",interface.gli_term_prev[subdomain]['wetting'])
...@@ -538,7 +542,8 @@ class DomainPatch(df.SubDomain): ...@@ -538,7 +542,8 @@ class DomainPatch(df.SubDomain):
# read neighbouring pressure values from interface and write them to # read neighbouring pressure values from interface and write them to
# the function self.neighbouring_interface_pressure # the function self.neighbouring_interface_pressure
interface.write_pressure_dofs(to_function = self.neighbouring_interface_pressure[phase], # pa_prev = df.interpolate(df.Constant(0), self.function_space[phase])
interface.write_pressure_dofs(to_function = pa_prev,#self.neighbouring_interface_pressure[phase], #
interface_dofs = interface_dofs[phase],# interface_dofs = interface_dofs[phase],#
dof_to_vert_map = self.dof2vertex[phase],# dof_to_vert_map = self.dof2vertex[phase],#
local_to_parent_vertex_map = self.parent_mesh_index,# local_to_parent_vertex_map = self.parent_mesh_index,#
...@@ -548,7 +553,7 @@ class DomainPatch(df.SubDomain): ...@@ -548,7 +553,7 @@ class DomainPatch(df.SubDomain):
# use the above to calculate the gli update with # use the above to calculate the gli update with
# in the paper p^{i-1}_{3-l} # in the paper p^{i-1}_{3-l}
pa_prev = self.neighbouring_interface_pressure[phase] # pa_prev = self.neighbouring_interface_pressure[phase]
phi_a = self.testfunction[phase] phi_a = self.testfunction[phase]
gli_pressure_part = df.assemble(-2*dt*Lambda[phase]*pa_prev*phi_a*ds) gli_pressure_part = df.assemble(-2*dt*Lambda[phase]*pa_prev*phi_a*ds)
gli_p_part = gli_pressure_part.get_local() gli_p_part = gli_pressure_part.get_local()
...@@ -634,44 +639,14 @@ class DomainPatch(df.SubDomain): ...@@ -634,44 +639,14 @@ class DomainPatch(df.SubDomain):
self.function_space = dict() self.function_space = dict()
self.trialfunction = dict() self.trialfunction = dict()
self.testfunction = dict() self.testfunction = dict()
self.neighbouring_interface_pressure = dict() # self.neighbouring_interface_pressure = dict()
for phase in self.has_phases: for phase in self.has_phases:
self.function_space.update({phase: df.FunctionSpace(self.mesh, 'P', 1)}) self.function_space.update({phase: df.FunctionSpace(self.mesh, 'P', 1)})
self.trialfunction.update({phase: df.TrialFunction(self.function_space[phase])}) self.trialfunction.update({phase: df.TrialFunction(self.function_space[phase])})
self.testfunction.update({phase: df.TestFunction(self.function_space[phase])}) self.testfunction.update({phase: df.TestFunction(self.function_space[phase])})
self.neighbouring_interface_pressure.update( # self.neighbouring_interface_pressure.update(
{phase: df.interpolate(df.Constant(0), self.function_space[phase])} # {phase: df.interpolate(df.Constant(0), self.function_space[phase])}
) # )
print(f"Python id of self.neighbouring_interface_pressure[{phase}]", id(self.neighbouring_interface_pressure[phase]))
# if self.isRichards:
# self.function_space = {'wetting': df.FunctionSpace(self.mesh, 'P', 1)}
# self.trialfunction = {
# 'wetting': df.TrialFunction(self.function_space['wetting'])
# }
# self.testfunction = {
# 'wetting': df.TestFunction(self.function_space['wetting'])
# }
# self.neighbouring_interface_pressure = {
# 'wetting' : df.interpolate(df.Constant(0), self.function_space['wetting'])
# }
# else:
# self.function_space = {#
# 'wetting' : df.FunctionSpace(self.mesh, 'P', 1),#
# 'nonwetting' : df.FunctionSpace(self.mesh, 'P', 1)#
# }
# self.trialfunction = {
# 'wetting' : df.TrialFunction(self.function_space['wetting']),
# 'nonwetting' : df.TrialFunction(self.function_space['nonwetting'])
# }
# self.testfunction = {
# 'wetting' : df.TestFunction(self.function_space['wetting']),
# 'nonwetting' : df.TestFunction(self.function_space['nonwetting'])
# }
# self.neighbouring_interface_pressure = {#
# 'wetting' : df.interpolate(df.Constant(0), self.function_space['wetting']),#
# 'nonwetting' : df.interpolate(df.Constant(0), self.function_space['nonwetting'])
# }
def _init_dof_maps(self) -> None: def _init_dof_maps(self) -> None:
""" calculate dof maps and the vertex to parent map. """ calculate dof maps and the vertex to parent map.
......
...@@ -85,13 +85,13 @@ isRichards = { ...@@ -85,13 +85,13 @@ isRichards = {
} }
##################### MESH ##################################################### ##################### MESH #####################################################
mesh_resolution = 10 mesh_resolution = 40
##################### TIME ##################################################### ##################### TIME #####################################################
timestep_size = 2*0.0001 timestep_size = 4*0.0001
number_of_timesteps = 30 number_of_timesteps = 5000
# 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_analise = 7 number_of_timesteps_to_analyse = 11
starttime = 0 starttime = 0
...@@ -115,8 +115,8 @@ L = {# ...@@ -115,8 +115,8 @@ L = {#
lambda_param = {# lambda_param = {#
# subdom_num : lambda parameter for the L-scheme # subdom_num : lambda parameter for the L-scheme
1 : {'wetting' :160},# 1 : {'wetting' :100},#
2 : {'wetting' :160} 2 : {'wetting' :100}
} }
## relative permeabilty functions on subdomain 1 ## relative permeabilty functions on subdomain 1
...@@ -147,9 +147,14 @@ relative_permeability = {# ...@@ -147,9 +147,14 @@ relative_permeability = {#
2: subdomain2_rel_perm 2: subdomain2_rel_perm
} }
def saturation(pressure, subdomain_index): # this function needs to be monotonically decreasing in the capillary_pressure.
# since in the richards case pc = -pw, this becomes as a function of pw a monotonically
# INCREASING function like in our Richards-Richards paper. However, since we unify
# the treatment in the code for Richards and two-phase, we need the same requierment
# for both cases, two-phase and Richards.
def saturation(capillary_pressure, subdomain_index):
# inverse capillary pressure-saturation-relationship # inverse capillary pressure-saturation-relationship
return df.conditional(pressure < 0, 1/((1 - pressure)**(1/(subdomain_index + 1))), 1) return df.conditional(capillary_pressure > 0, 1/((1 + capillary_pressure)**(1/(subdomain_index + 1))), 1)
sat_pressure_relationship = {# sat_pressure_relationship = {#
1: ft.partial(saturation, subdomain_index = 1),# 1: ft.partial(saturation, subdomain_index = 1),#
...@@ -197,7 +202,7 @@ write_to_file = { ...@@ -197,7 +202,7 @@ write_to_file = {
} }
# initialise LDD simulation class # initialise LDD simulation class
simulation = ldd.LDDsimulation(tol = 1E-12) simulation = ldd.LDDsimulation(tol = 1E-14)
simulation.set_parameters(output_dir = "./output/",# simulation.set_parameters(output_dir = "./output/",#
subdomain_def_points = subdomain_def_points,# subdomain_def_points = subdomain_def_points,#
isRichards = isRichards,# isRichards = isRichards,#
...@@ -213,11 +218,11 @@ simulation.set_parameters(output_dir = "./output/",# ...@@ -213,11 +218,11 @@ simulation.set_parameters(output_dir = "./output/",#
saturation = sat_pressure_relationship,# saturation = sat_pressure_relationship,#
starttime = starttime,# starttime = starttime,#
number_of_timesteps = number_of_timesteps, number_of_timesteps = number_of_timesteps,
number_of_timesteps_to_analise = number_of_timesteps_to_analise, number_of_timesteps_to_analyse = number_of_timesteps_to_analyse,
timestep_size = timestep_size,# timestep_size = timestep_size,#
sources = source_expression,# sources = source_expression,#
initial_conditions = initial_condition,# initial_conditions = initial_condition,#
dirichletBC_Expressions = dirichletBC,# dirichletBC_expression_strings = dirichletBC,#
exact_solution = exact_solution,# exact_solution = exact_solution,#
write2file = write_to_file,# write2file = write_to_file,#
) )
......
...@@ -89,9 +89,9 @@ isRichards = { ...@@ -89,9 +89,9 @@ isRichards = {
############ GRID ########################ü ############ GRID ########################ü
mesh_resolution = 100 mesh_resolution = 40
timestep_size = 2*0.001 timestep_size = 0.001
number_of_timesteps = 1000 number_of_timesteps = 2000
# 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
...@@ -122,9 +122,9 @@ L = {# ...@@ -122,9 +122,9 @@ L = {#
lambda_param = {# lambda_param = {#
# subdom_num : lambda parameter for the L-scheme # subdom_num : lambda parameter for the L-scheme
1 : {'wetting' :140, 1 : {'wetting' :140,
'nonwetting': 1000},# 'nonwetting': 2400},#
2 : {'wetting' :140, 2 : {'wetting' :140,
'nonwetting': 1000} 'nonwetting': 2400}
} }
## relative permeabilty functions on subdomain 1 ## relative permeabilty functions on subdomain 1
...@@ -147,7 +147,7 @@ def rel_perm2w(s): ...@@ -147,7 +147,7 @@ def rel_perm2w(s):
# relative permeabilty wetting on subdomain2 # relative permeabilty wetting on subdomain2
return s**3 return s**3
def rel_perm2nw(s): def rel_perm2nw(s):
# relative permeabilty nonwetting on subdomain2 # relative permeabilty nonwetting on subdosym.cos(0.8*t - (0.8*x + 1/7*y))main2
return (1-s)**2 return (1-s)**2
_rel_perm2w = ft.partial(rel_perm2w) _rel_perm2w = ft.partial(rel_perm2w)
...@@ -168,7 +168,7 @@ relative_permeability = {# ...@@ -168,7 +168,7 @@ relative_permeability = {#
# we set alpha = 0, assume m = 1-1/n (see Helmig) and assume that residual saturation is Sw # we set alpha = 0, assume m = 1-1/n (see Helmig) and assume that residual saturation is Sw
def saturation(capillary_pressure, n_index, alpha): def saturation(capillary_pressure, n_index, alpha):
# inverse capillary pressure-saturation-relationship # inverse capillary pressure-saturation-relationship
return df.conditional(capillary_pressure < 0, 1/((1 + (alpha*capillary_pressure)**n_index)**((n_index - 1)/n_index)), 1) return df.conditional(capillary_pressure > 0, 1/((1 + (alpha*capillary_pressure)**n_index)**((n_index - 1)/n_index)), 1)
# S-pc-relation ship. We use the van Genuchten approach, i.e. pc = 1/alpha*(S^{-1/m} -1)^1/n, where # S-pc-relation ship. We use the van Genuchten approach, i.e. pc = 1/alpha*(S^{-1/m} -1)^1/n, where
# we set alpha = 0, assume m = 1-1/n (see Helmig) and assume that residual saturation is Sw # we set alpha = 0, assume m = 1-1/n (see Helmig) and assume that residual saturation is Sw
...@@ -318,7 +318,7 @@ write_to_file = { ...@@ -318,7 +318,7 @@ write_to_file = {
# initialise LDD simulation class # initialise LDD simulation class
simulation = ldd.LDDsimulation(tol = 1E-12) simulation = ldd.LDDsimulation(tol = 1E-14)
simulation.set_parameters(output_dir = "./output/",# simulation.set_parameters(output_dir = "./output/",#
subdomain_def_points = subdomain_def_points,# subdomain_def_points = subdomain_def_points,#
isRichards = isRichards,# isRichards = isRichards,#
...@@ -338,7 +338,7 @@ simulation.set_parameters(output_dir = "./output/",# ...@@ -338,7 +338,7 @@ simulation.set_parameters(output_dir = "./output/",#
timestep_size = timestep_size,# timestep_size = timestep_size,#
sources = source_expression,# sources = source_expression,#
initial_conditions = initial_condition,# initial_conditions = initial_condition,#
dirichletBC_Expressions = dirichletBC,# dirichletBC_expression_strings = dirichletBC,#
exact_solution = exact_solution,# exact_solution = exact_solution,#
write2file = write_to_file,# write2file = write_to_file,#
) )
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment