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

update RR-multi-patch-with-inner-patch example and run as test

parent 23cb8e09
No related branches found
No related tags found
No related merge requests found
......@@ -7,28 +7,34 @@ import sympy as sym
# import domainPatch as dp
import LDDsimulation as ldd
import functools as ft
import helpers as hlp
# import ufl as ufl
# init sympy session
sym.init_printing()
# ----------------------------------------------------------------------------#
# ------------------- MESH ---------------------------------------------------#
# ----------------------------------------------------------------------------#
mesh_resolution = 51
# ----------------------------------------:-------------------------------------#
# ------------------- TIME ---------------------------------------------------#
# ----------------------------------------------------------------------------#
use_case = "RR-multi-patch-with-inner-patch-cuttoff"
solver_tol = 1E-6
############ GRID #######################ü
mesh_resolution = 20
timestep_size = 0.01
number_of_timesteps = 150
# decide how many timesteps you want analysed. Analysed means, that we write
# out subsequent errors of the L-iteration within the timestep.
# decide how many timesteps you want analysed. Analysed means, that we write out
# subsequent errors of the L-iteration within the timestep.
number_of_timesteps_to_analyse = 10
starttime = 0
# Lw = 0.25/timestep_size
Lw = 0.25
Lw = 0.25 #/timestep_size
Lnw=Lw
lambda_w = 4
include_gravity = True
debugflag = False
analyse_condition = False
output_string = "./output/new_gravity_cutoff_term-number_of_timesteps{}_".format(number_of_timesteps)
# ----------------------------------------------------------------------------#
# ------------------- Domain and Interface -----------------------------------#
......@@ -205,15 +211,15 @@ L = {
5: {'wetting': Lw}
}
lamdal_w = 32
lambda_w = 32
# subdom_num : lambda parameter for the L-scheme
lambda_param = {
1: {'wetting': lamdal_w},
2: {'wetting': lamdal_w},
3: {'wetting': lamdal_w},
4: {'wetting': lamdal_w},
5: {'wetting': lamdal_w}
1: {'wetting': lambda_w},
2: {'wetting': lambda_w},
3: {'wetting': lambda_w},
4: {'wetting': lambda_w},
5: {'wetting': lambda_w}
}
......@@ -460,84 +466,37 @@ p_e_sym = {
}
pc_e_sym = {
1: -1*p_e_sym[1]['wetting'],
2: -1*p_e_sym[2]['wetting'],
3: -1*p_e_sym[3]['wetting'],
4: -1*p_e_sym[4]['wetting'],
5: -1*p_e_sym[5]['wetting']
}
# turn above symbolic code into exact solution for dolphin and
# construct the rhs that matches the above exact solution.
dtS = dict()
div_flux = dict()
source_expression = dict()
exact_solution = dict()
initial_condition = dict()
pc_e_sym = dict()
for subdomain, isR in isRichards.items():
dtS.update({subdomain: dict()})
div_flux.update({subdomain: dict()})
source_expression.update({subdomain: dict()})
exact_solution.update({subdomain: dict()})
initial_condition.update({subdomain: dict()})
if isR:
subdomain_has_phases = ["wetting"]
pc_e_sym.update({subdomain: -p_e_sym[subdomain]['wetting']})
else:
subdomain_has_phases = ["wetting", "nonwetting"]
# conditional for S_pc_prime
pc = pc_e_sym[subdomain]
dtpc = sym.diff(pc, t, 1)
dxpc = sym.diff(pc, x, 1)
dypc = sym.diff(pc, y, 1)
S = sym.Piecewise((S_pc_sym[subdomain](pc), pc > 0), (1, True))
dS = sym.Piecewise((S_pc_sym_prime[subdomain](pc), pc > 0), (0, True))
for phase in subdomain_has_phases:
# Turn above symbolic expression for exact solution into c code
exact_solution[subdomain].update(
{phase: sym.printing.ccode(p_e_sym[subdomain][phase])}
)
# save the c code for initial conditions
initial_condition[subdomain].update(
{phase: sym.printing.ccode(p_e_sym[subdomain][phase].subs(t, 0))}
)
if phase == "nonwetting":
dtS[subdomain].update(
{phase: -porosity[subdomain]*dS*dtpc}
)
else:
dtS[subdomain].update(
{phase: porosity[subdomain]*dS*dtpc}
)
pa = p_e_sym[subdomain][phase]
dxpa = sym.diff(pa, x, 1)
dxdxpa = sym.diff(pa, x, 2)
dypa = sym.diff(pa, y, 1)
dydypa = sym.diff(pa, y, 2)
mu = viscosity[subdomain][phase]
ka = relative_permeability[subdomain][phase]
dka = ka_prime[subdomain][phase]
rho = densities[subdomain][phase]
g = gravity_acceleration
if phase == "nonwetting":
# x part of div(flux) for nonwetting
dxdxflux = -1/mu*dka(1-S)*dS*dxpc*dxpa + 1/mu*dxdxpa*ka(1-S)
# y part of div(flux) for nonwetting
dydyflux = -1/mu*dka(1-S)*dS*dypc*(dypa - rho*g) \
+ 1/mu*dydypa*ka(1-S)
else:
# x part of div(flux) for wetting
dxdxflux = 1/mu*dka(S)*dS*dxpc*dxpa + 1/mu*dxdxpa*ka(S)
# y part of div(flux) for wetting
dydyflux = 1/mu*dka(S)*dS*dypc*(dypa - rho*g) + 1/mu*dydypa*ka(S)
div_flux[subdomain].update({phase: dxdxflux + dydyflux})
contructed_rhs = dtS[subdomain][phase] - div_flux[subdomain][phase]
source_expression[subdomain].update(
{phase: sym.printing.ccode(contructed_rhs)}
pc_e_sym.update({subdomain: p_e_sym[subdomain]['nonwetting']
- p_e_sym[subdomain]['wetting']})
symbols = {"x": x,
"y": y,
"t": t}
# turn above symbolic code into exact solution for dolphin and
# construct the rhs that matches the above exact solution.
exact_solution_example = hlp.generate_exact_solution_expressions(
symbols=symbols,
isRichards=isRichards,
symbolic_pressure=p_e_sym,
symbolic_capillary_pressure=pc_e_sym,
saturation_pressure_relationship=S_pc_sym,
saturation_pressure_relationship_prime=S_pc_sym_prime,
viscosity=viscosity,
porosity=porosity,
relative_permeability=relative_permeability,
relative_permeability_prime=ka_prime,
densities=densities,
gravity_acceleration=gravity_acceleration,
include_gravity=include_gravity,
)
# print(f"source_expression[{subdomain}][{phase}] =", source_expression[subdomain][phase])
source_expression = exact_solution_example['source']
exact_solution = exact_solution_example['exact_solution']
initial_condition = exact_solution_example['initial_condition']
# Dictionary of dirichlet boundary conditions.
dirichletBC = dict()
......@@ -568,87 +527,15 @@ for subdomain in isRichards.keys():
)
# # construction of the rhs that matches the above exact solution.
# dtS = dict()
# div_flux = dict()
# source_expression = dict()
# for subdomain, isR in isRichards.items():
# dtS.update({subdomain: dict()})
# div_flux.update({subdomain: dict()})
# source_expression.update({subdomain: dict()})
# if isR:
# subdomain_has_phases = ["wetting"]
# else:
# subdomain_has_phases = ["wetting", "nonwetting"]
#
# # conditional for S_pc_prime
# pc = pc_e_sym[subdomain]
# dtpc = sym.diff(pc, t, 1)
# dxpc = sym.diff(pc, x, 1)
# dypc = sym.diff(pc, y, 1)
# S = sym.Piecewise((S_pc_sym[subdomain](pc), pc > 0), (1, True))
# dS = sym.Piecewise((S_pc_sym_prime[subdomain](pc), pc > 0), (0, True))
# for phase in subdomain_has_phases:
# if phase == "nonwetting":
# dS = -dS
# dtS[subdomain].update(
# {phase: porosity[subdomain]*dS*dtpc}
# )
# pa = p_e_sym[subdomain][phase]
# dxpa = sym.diff(pa, x, 1)
# dxdxpa = sym.diff(pa, x, 2)
# dypa = sym.diff(pa, y, 1)
# dydypa = sym.diff(pa, y, 2)
# mu = viscosity[subdomain][phase]
# ka = relative_permeability[subdomain][phase]
# dka = ka_prime[subdomain][phase]
# rho = densities[subdomain][phase]
# g = gravity_acceleration
#
# if phase == "nonwetting":
# dxdxflux = 1/mu*dka(1-S)*dS*dxpc*dxpa - 1/mu*dxdxpa*ka(1-S)
# dydyflux = 1/mu*dka(1-S)*dS*dypc*(dypa - rho*g) \
# - 1/mu*dydypa*ka(1-S)
# else:
# dxdxflux = -1/mu*dka(S)*dS*dxpc*dxpa - 1/mu*dxdxpa*ka(S)
# dydyflux = -1/mu*dka(S)*dS*dypc*(dypa - rho*g) - 1/mu*dydypa*ka(S)
# div_flux[subdomain].update({phase: dxdxflux + dydyflux})
# contructed_rhs = dtS[subdomain][phase] + div_flux[subdomain][phase]
# source_expression[subdomain].update(
# {phase: sym.printing.ccode(contructed_rhs)}
# )
# print(f"source_expression[{subdomain}][{phase}] =",
# source_expression[subdomain][phase])
# similarly to the outer boundary dictionary, if a patch has no outer boundary
# None should be written instead of an expression.
# This is a bit of a brainfuck:
# dirichletBC[ind] gives a dictionary of the outer boundaries of subdomain ind.
# Since a domain patch can have several disjoint outer boundary parts, the
# expressions need to get an enumaration index which starts at 0.
# So dirichletBC[ind][j] is the dictionary of outer dirichlet conditions of
# subdomain ind and boundary part j.
# Finally, dirichletBC[ind][j]['wetting'] and dirichletBC[ind][j]['nonwetting']
# return the actual expression needed for the dirichlet condition for both
# phases if present.
# subdomain index: {outer boudary part index: {phase: expression}}
# dirichletBC = {
# 1: {0: {'wetting': exact_solution[1]['wetting']}},
# 2: {0: {'wetting': exact_solution[2]['wetting']}},
# 3: {0: {'wetting': exact_solution[3]['wetting']}},
# 4: {0: {'wetting': exact_solution[4]['wetting']}},
# 5: {0: {'wetting': exact_solution[5]['wetting']}}
# }
write_to_file = {
'meshes_and_markers': True,
'L_iterations': True
}
# initialise LDD simulation class
simulation = ldd.LDDsimulation(tol=1E-14, debug=False, LDDsolver_tol=1E-7)
simulation.set_parameters(output_dir="./output/with_cutoff_function",
simulation = ldd.LDDsimulation(tol=1E-14, debug=debugflag, LDDsolver_tol=solver_tol)
simulation.set_parameters(use_case=use_case,
output_dir=output_string,
subdomain_def_points=subdomain_def_points,
isRichards=isRichards,
interface_def_points=interface_def_points,
......@@ -670,12 +557,12 @@ simulation.set_parameters(output_dir="./output/with_cutoff_function",
dirichletBC_expression_strings=dirichletBC,
exact_solution=exact_solution,
densities=densities,
include_gravity=True,
include_gravity=include_gravity,
write2file=write_to_file,
)
simulation.initialise()
# print(simulation.__dict__)
simulation.run()
simulation.run(analyse_condition=analyse_condition)
# simulation.LDDsolver(time=0, debug=True, analyse_timestep=True)
# df.info(parameters, True)
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment