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

setup layered soil simulations

parent 9ffa5580
Branches
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 1001 additions and 127 deletions
Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-all-params-one.py
+ 3
− 2
View file @ d3ff68d5
......@@ -24,6 +24,8 @@ datestr = date.strftime("%Y-%m-%d")
sym.init_printing()
# solver_tol = 6E-7
use_case = "TP-R-layered-soil-all-params-set-one-new-lambda"
thisfile = "TP-R-layered_soil-all-params-one.py"
max_iter_num = 500
FEM_Lagrange_degree = 1
mesh_study = False
......@@ -603,11 +605,10 @@ for subdomain in isRichards.keys():
# read this file and print it to std out. This way the simulation can produce a
# log file with ./TP-R-layered_soil.py | tee simulation.log
f = open('TP-R-layered_soil.py', 'r')
f = open(thisfile, 'r')
print(f.read())
f.close()
for starttime in starttimes:
for mesh_resolution, solver_tol in resolutions.items():
# initialise LDD simulation class
......
Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil.py
+ 44
− 31
View file @ d3ff68d5
......@@ -23,32 +23,34 @@ datestr = date.strftime("%Y-%m-%d")
# init sympy session
sym.init_printing()
# solver_tol = 6E-7
use_case = "TP-R-layered-soil-realistic-new-lambda"
max_iter_num = 600
use_case = "TP-R-layered-soil-realistic"
# name of this very file. Needed for log output.
thisfile = "TP-R-layered_soil.py"
max_iter_num = 300
FEM_Lagrange_degree = 1
mesh_study = False
resolutions = {
# 1: 1e-7, # h=2
# 2: 2e-5, # h=1.1180
# 4: 1e-6, # h=0.5590
# 8: 1e-6, # h=0.2814
# 16: 1e-6, # h=0.1412
32: 1e-6,
# 64: 5e-7,
# 128: 5e-7
# 1: 2e-6, # h=2
# 2: 2e-6, # h=1.1180
# 4: 2e-6, # h=0.5590
# 8: 2e-6, # h=0.2814
# 16: 2e-6, # h=0.1412
32: 2e-6,
# 64: 2e-6,
# 128: 2e-6
}
# GRID #######################
# mesh_resolution = 20
timestep_size = 0.005
number_of_timesteps = 200
plot_timestep_every = 2
timestep_size = 0.001
number_of_timesteps = 1000
plot_timestep_every = 4
# 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 = 5
number_of_timesteps_to_analyse = 4
starttimes = [0.0]
Lw1 = 0.25 # /timestep_size
Lw1 = 0.025 # /timestep_size
Lnw1 = Lw1
Lw2 = 0.025 # /timestep_size
Lnw2 = Lw2
......@@ -57,8 +59,8 @@ Lnw3 = Lw3
Lw4 = 0.025 # /timestep_size
Lnw4 = Lw4
lambda12_w = 4
lambda12_nw = 4
lambda12_w = 40
lambda12_nw = 40
lambda23_w = 40
lambda23_nw = 40
lambda34_w = 40
......@@ -66,7 +68,7 @@ lambda34_nw = 40
include_gravity = False
debugflag = False
analyse_condition = False
analyse_condition = True
if mesh_study:
output_string = "./output/{}-{}_timesteps{}_P{}".format(
......@@ -279,18 +281,28 @@ viscosity = {
}
# Dict of the form: { subdom_num : density }
# densities = {
# 1: {'wetting': 9.97, # 997
# 'nonwetting': 0.01225}, # 1.225}},
# 2: {'wetting': 9.97, # 997
# 'nonwetting': 0.01225}, # 1.225}},
# 3: {'wetting': 9.97, # 997
# 'nonwetting': 0.01225}, # 1.225}},
# 4: {'wetting': 9.97, # 997
# 'nonwetting': 0.01225}, # 1.225}}
# }
densities = {
1: {'wetting': 9.97, # 997
'nonwetting': 0.01225}, # 1.225}},
2: {'wetting': 9.97, # 997
'nonwetting': 0.01225}, # 1.225}},
3: {'wetting': 9.97, # 997
'nonwetting': 0.01225}, # 1.225}},
4: {'wetting': 9.97, # 997
'nonwetting': 0.01225}, # 1.225}}
1: {'wetting': 997, # 997
'nonwetting': 1.225}, # 1.225}},
2: {'wetting': 997, # 997
'nonwetting': 1.225}, # 1.225}},
3: {'wetting': 997, # 997
'nonwetting': 1.225}, # 1.225}},
4: {'wetting': 997, # 997
'nonwetting': 1.225}, # 1.225}}
}
gravity_acceleration = 9.81
# porosities taken from
# https://www.geotechdata.info/parameter/soil-porosity.html
......@@ -338,10 +350,10 @@ lambda_param = {
# after Lewis, see pdf file
intrinsic_permeability = {
1: 1, # sand
1: 0.1, # sand
2: 0.1, # sand, there is a range
3: 10e-2, # clay has a range
4: 10e-4
3: 0.001, #10e-2, # clay has a range
4: 0.001, #10e-3
}
......@@ -698,9 +710,10 @@ for subdomain in isRichards.keys():
{outer_boundary_ind: exact_solution[subdomain]}
)
# read this file and print it to std out. This way the simulation can produce a
# log file with ./TP-R-layered_soil.py | tee simulation.log
f = open('TP-R-layered_soil.py', 'r')
f = open(thisfile, 'r')
print(f.read())
f.close()
......
Two-phase-Richards/multi-patch/layered_soil/mesh_study/TP-R-layered_soil.py 0 → 100755
+ 803
− 0
View file @ d3ff68d5
This diff is collapsed.
Two-phase-Richards/multi-patch/layered_soil/mesh_study/run-simulation 0 → 100755
+ 16
− 0
View file @ d3ff68d5
#!/bin/bash
[ $# -eq 0 ] && { echo "Usage: $0 simulation_file [logfile_name]"; exit 1; }
SIMULATION_FILE=$1
SIMULATION=${SIMULATION_FILE%.py}
LOGFILE_DEFAULT="$SIMULATION.log"
DATE=$(date -I)
LOGFILE=${2:-$DATE-$LOGFILE_DEFAULT}
GREETING="Simulation $SIMULATION is run on $DATE by $USER"
echo $GREETING
echo "running $SIMULATION_FILE | tee $LOGFILE"
./$SIMULATION_FILE | tee $LOGFILE
Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-test.py
+ 70
− 41
View file @ d3ff68d5
......@@ -20,6 +20,7 @@ datestr = date.strftime("%Y-%m-%d")
sym.init_printing()
use_case = "TP-R-2-patch-realistic-new-implementation"
thisfile = "TP-R-2-patch-test.py"
# solver_tol = 6E-7
max_iter_num = 1000
FEM_Lagrange_degree = 1
......@@ -56,7 +57,17 @@ include_gravity = True
debugflag = True
analyse_condition = True
output_string = "./output/{}-{}_timesteps{}_P{}".format(datestr, use_case, number_of_timesteps, FEM_Lagrange_degree)
if mesh_study:
output_string = "./output/{}-{}_timesteps{}_P{}".format(
datestr, use_case, number_of_timesteps, FEM_Lagrange_degree
)
else:
for tol in resolutions.values():
solver_tol = tol
output_string = "./output/{}-{}_timesteps{}_P{}_solver_tol{}".format(
datestr, use_case, number_of_timesteps, FEM_Lagrange_degree, solver_tol
)
# toggle what should be written to files
if mesh_study:
......@@ -433,6 +444,10 @@ for subdomain in isRichards.keys():
{outer_boundary_ind: exact_solution[subdomain]}
)
f = open(thisfile, 'r')
print(f.read())
f.close()
for starttime in starttimes:
for mesh_resolution, solver_tol in resolutions.items():
......@@ -446,33 +461,34 @@ for starttime in starttimes:
mesh_study=mesh_study
)
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,
outer_boundary_def_points=outer_boundary_def_points,
adjacent_subdomains=adjacent_subdomains,
mesh_resolution=mesh_resolution,
viscosity=viscosity,
porosity=porosity,
L=L,
lambda_param=lambda_param,
relative_permeability=relative_permeability,
saturation=sat_pressure_relationship,
starttime=starttime,
number_of_timesteps=number_of_timesteps,
number_of_timesteps_to_analyse=number_of_timesteps_to_analyse,
plot_timestep_every=plot_timestep_every,
timestep_size=timestep_size,
sources=source_expression,
initial_conditions=initial_condition,
dirichletBC_expression_strings=dirichletBC,
exact_solution=exact_solution,
densities=densities,
include_gravity=include_gravity,
write2file=write_to_file,
)
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,
outer_boundary_def_points=outer_boundary_def_points,
adjacent_subdomains=adjacent_subdomains,
mesh_resolution=mesh_resolution,
viscosity=viscosity,
porosity=porosity,
L=L,
lambda_param=lambda_param,
relative_permeability=relative_permeability,
saturation=sat_pressure_relationship,
starttime=starttime,
number_of_timesteps=number_of_timesteps,
number_of_timesteps_to_analyse=number_of_timesteps_to_analyse,
plot_timestep_every=plot_timestep_every,
timestep_size=timestep_size,
sources=source_expression,
initial_conditions=initial_condition,
dirichletBC_expression_strings=dirichletBC,
exact_solution=exact_solution,
densities=densities,
include_gravity=include_gravity,
write2file=write_to_file,
)
simulation.initialise()
output_dir = simulation.output_dir
......@@ -480,26 +496,39 @@ for starttime in starttimes:
output = simulation.run(analyse_condition=analyse_condition)
for subdomain_index, subdomain_output in output.items():
mesh_h = subdomain_output['mesh_size']
for phase, different_errornorms in subdomain_output['errornorm'].items():
filename = output_dir + "subdomain{}-space-time-errornorm-{}-phase.csv".format(subdomain_index, phase)
# for errortype, errornorm in different_errornorms.items():
# eocfile = open("eoc_filename", "a")
# eocfile.write( str(mesh_h) + " " + str(errornorm) + "\n" )
# eocfile.close()
# if subdomain.isRichards:mesh_h
for phase, error_dict in subdomain_output['errornorm'].items():
filename = output_dir \
+ "subdomain{}".format(subdomain_index)\
+ "-space-time-errornorm-{}-phase.csv".format(phase)
# for errortype, errornorm in error_dict.items():
# eocfile = open("eoc_filename", "a")
# eocfile.write( str(mesh_h) + " " + str(errornorm) + "\n" )
# eocfile.close()
# if subdomain.isRichards:mesh_h
data_dict = {
'mesh_parameter': mesh_resolution,
'mesh_h': mesh_h,
}
for error_type, errornorms in different_errornorms.items():
for norm_type, errornorm in error_dict.items():
data_dict.update(
{error_type: errornorms}
{norm_type: errornorm}
)
errors = pd.DataFrame(data_dict, index=[mesh_resolution])
# check if file exists
if os.path.isfile(filename) == True:
if os.path.isfile(filename) is True:
with open(filename, 'a') as f:
errors.to_csv(f, header=False, sep='\t', encoding='utf-8', index=False)
errors.to_csv(
f,
header=False,
sep='\t',
encoding='utf-8',
index=False
)
else:
errors.to_csv(filename, sep='\t', encoding='utf-8', index=False)
errors.to_csv(
filename,
sep='\t',
encoding='utf-8',
index=False
)
Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-realistic-parameters-densities-scaled.py
+ 65
− 53
View file @ d3ff68d5
......@@ -20,6 +20,7 @@ datestr = date.strftime("%Y-%m-%d")
sym.init_printing()
use_case = "TPR-2-desities-scaled-down"
thisfile = "TP-R-2-realistic-parameters-densities-scaled.py"
# solver_tol = 6E-7
max_iter_num = 500
FEM_Lagrange_degree = 1
......@@ -57,11 +58,15 @@ debugflag = False
analyse_condition = True
if mesh_study:
output_string = "./output/{}-{}_timesteps{}_P{}".format(datestr, use_case, number_of_timesteps, FEM_Lagrange_degree)
output_string = "./output/{}-{}_timesteps{}_P{}".format(
datestr, use_case, number_of_timesteps, FEM_Lagrange_degree
)
else:
for tol in resolutions.values():
solver_tol = tol
output_string = "./output/{}-{}_timesteps{}_P{}_solver_tol{}".format(datestr, use_case, number_of_timesteps, FEM_Lagrange_degree, solver_tol)
output_string = "./output/{}-{}_timesteps{}_P{}_solver_tol{}".format(
datestr, use_case, number_of_timesteps, FEM_Lagrange_degree, solver_tol
)
# toggle what should be written to files
if mesh_study:
......@@ -188,19 +193,12 @@ L = {#
}
lambda_param = {#
# subdom_num : lambda parameter for the L-scheme
lambda_param = {
# subdom_num : lambda parameter for the L-scheme
0 : {'wetting' :lambda_w,
'nonwetting': lambda_nw},#
'nonwetting': lambda_nw},
}
# intrinsic_permeability = {
# 1: {"wetting": 1,
# "nonwetting": 1},
# 2: {"wetting": 1,
# "nonwetting": 1},
# }
intrinsic_permeability = {
1: 1,
2: 1,
......@@ -459,7 +457,7 @@ for subdomain in isRichards.keys():
#
# sa
f = open('TP-R-2-patch-mesh-study.py', 'r')
f = open(thisfile, 'r')
print(f.read())
f.close()
......@@ -476,33 +474,34 @@ for starttime in starttimes:
mesh_study=mesh_study
)
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,
outer_boundary_def_points=outer_boundary_def_points,
adjacent_subdomains=adjacent_subdomains,
mesh_resolution=mesh_resolution,
viscosity=viscosity,
porosity=porosity,
L=L,
lambda_param=lambda_param,
relative_permeability=relative_permeability,
saturation=sat_pressure_relationship,
starttime=starttime,
number_of_timesteps=number_of_timesteps,
number_of_timesteps_to_analyse=number_of_timesteps_to_analyse,
plot_timestep_every=plot_timestep_every,
timestep_size=timestep_size,
sources=source_expression,
initial_conditions=initial_condition,
dirichletBC_expression_strings=dirichletBC,
exact_solution=exact_solution,
densities=densities,
include_gravity=include_gravity,
write2file=write_to_file,
)
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,
outer_boundary_def_points=outer_boundary_def_points,
adjacent_subdomains=adjacent_subdomains,
mesh_resolution=mesh_resolution,
viscosity=viscosity,
porosity=porosity,
L=L,
lambda_param=lambda_param,
relative_permeability=relative_permeability,
saturation=sat_pressure_relationship,
starttime=starttime,
number_of_timesteps=number_of_timesteps,
number_of_timesteps_to_analyse=number_of_timesteps_to_analyse,
plot_timestep_every=plot_timestep_every,
timestep_size=timestep_size,
sources=source_expression,
initial_conditions=initial_condition,
dirichletBC_expression_strings=dirichletBC,
exact_solution=exact_solution,
densities=densities,
include_gravity=include_gravity,
write2file=write_to_file,
)
simulation.initialise()
output_dir = simulation.output_dir
......@@ -510,26 +509,39 @@ for starttime in starttimes:
output = simulation.run(analyse_condition=analyse_condition)
for subdomain_index, subdomain_output in output.items():
mesh_h = subdomain_output['mesh_size']
for phase, different_errornorms in subdomain_output['errornorm'].items():
filename = output_dir + "subdomain{}-space-time-errornorm-{}-phase.csv".format(subdomain_index, phase)
# for errortype, errornorm in different_errornorms.items():
# eocfile = open("eoc_filename", "a")
# eocfile.write( str(mesh_h) + " " + str(errornorm) + "\n" )
# eocfile.close()
# if subdomain.isRichards:mesh_h
for phase, error_dict in subdomain_output['errornorm'].items():
filename = output_dir \
+ "subdomain{}".format(subdomain_index)\
+ "-space-time-errornorm-{}-phase.csv".format(phase)
# for errortype, errornorm in error_dict.items():
# eocfile = open("eoc_filename", "a")
# eocfile.write( str(mesh_h) + " " + str(errornorm) + "\n" )
# eocfile.close()
# if subdomain.isRichards:mesh_h
data_dict = {
'mesh_parameter': mesh_resolution,
'mesh_h': mesh_h,
}
for error_type, errornorms in different_errornorms.items():
for norm_type, errornorm in error_dict.items():
data_dict.update(
{error_type: errornorms}
{norm_type: errornorm}
)
errors = pd.DataFrame(data_dict, index=[mesh_resolution])
# check if file exists
if os.path.isfile(filename) == True:
if os.path.isfile(filename) is True:
with open(filename, 'a') as f:
errors.to_csv(f, header=False, sep='\t', encoding='utf-8', index=False)
errors.to_csv(
f,
header=False,
sep='\t',
encoding='utf-8',
index=False
)
else:
errors.to_csv(filename, sep='\t', encoding='utf-8', index=False)
errors.to_csv(
filename,
sep='\t',
encoding='utf-8',
index=False
)
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment