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

update and clean up example

parent 19afee6c
Branches
No related tags found
No related merge requests found
#!/usr/bin/python3 #!/usr/bin/python3
"""TPR 2 patch soil simulation.
This program sets up an LDD simulation
"""
import dolfin as df import dolfin as df
import mshr
import numpy as np
import sympy as sym import sympy as sym
import typing as tp
import domainPatch as dp
import LDDsimulation as ldd
import functools as ft import functools as ft
import LDDsimulation as ldd
import helpers as hlp import helpers as hlp
import datetime import datetime
import os import os
import pandas as pd import pandas as pd
date = datetime.datetime.now()
datestr = date.strftime("%Y-%m-%d")
#import ufl as ufl
# init sympy session # init sympy session
sym.init_printing() sym.init_printing()
# PREREQUISITS ###############################################################
# check if output directory "./output" exists. This will be used in
# the generation of the output string.
if not os.path.exists('./output'):
os.mkdir('./output')
print("Directory ", './output', " created ")
else:
print("Directory ", './output', " already exists. Will use as output \
directory")
date = datetime.datetime.now()
datestr = date.strftime("%Y-%m-%d")
# Name of the usecase that will be printed during simulation.
use_case = "TPR-2-patch-realistic-testrun" use_case = "TPR-2-patch-realistic-testrun"
# solver_tol = 6E-7 # The name of this very file. Needed for creating log output.
thisfile = "TP-R-2-patch-mesh-study.py"
# GENERAL SOLVER CONFIG ######################################################
# maximal iteration per timestep
max_iter_num = 500 max_iter_num = 500
FEM_Lagrange_degree = 1 FEM_Lagrange_degree = 1
# GRID AND MESH STUDY SPECIFICATIONS #########################################
mesh_study = True mesh_study = True
resolutions = { resolutions = {
1: 1e-6, 1: 1e-6,
...@@ -36,18 +53,19 @@ resolutions = { ...@@ -36,18 +53,19 @@ resolutions = {
256: 1e-6, 256: 1e-6,
} }
############ GRID ####################### # starttimes gives a list of starttimes to run the simulation from.
# mesh_resolution = 20 # The list is looped over and a simulation is run with t_0 as initial time
# for each element t_0 in starttimes.
starttimes = [0.0]
timestep_size = 0.001 timestep_size = 0.001
number_of_timesteps = 1000 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
starttimes = [0.0]
Lw = 0.025 #/timestep_size # LDD scheme parameters ######################################################
Lnw= 0.025 Lw1 = 0.025 #/timestep_size
Lnw1= 0.025
Lw2 = 0.025 #/timestep_size
Lnw2= 0.025
lambda_w = 40 lambda_w = 40
lambda_nw = 40 lambda_nw = 40
...@@ -56,22 +74,38 @@ include_gravity = False ...@@ -56,22 +74,38 @@ include_gravity = False
debugflag = False debugflag = False
analyse_condition = True analyse_condition = True
if mesh_study: # I/O CONFIG #################################################################
output_string = "./output/{}-{}_timesteps{}_P{}".format(datestr, use_case, number_of_timesteps, FEM_Lagrange_degree) # when number_of_timesteps is high, it might take a long time to write all
else: # timesteps to disk. Therefore, you can choose to only write data of every
for tol in resolutions.values(): # plot_timestep_every timestep to disk.
solver_tol = tol plot_timestep_every = 4
output_string = "./output/{}-{}_timesteps{}_P{}_solver_tol{}".format(datestr, use_case, number_of_timesteps, FEM_Lagrange_degree, solver_tol) # Decide how many timesteps you want analysed. Analysed means, that
# subsequent errors of the L-iteration within the timestep are written out.
number_of_timesteps_to_analyse = 5
# toggle what should be written to files # fine grained control over data to be written to disk in the mesh study case
# as well as for a regular simuation for a fixed grid.
if mesh_study: if mesh_study:
write_to_file = { write_to_file = {
# output the relative errornorm (integration in space) w.r.t. an exact
# solution for each timestep into a csv file.
'space_errornorms': True, 'space_errornorms': True,
# save the mesh and marker functions to disk
'meshes_and_markers': True, 'meshes_and_markers': True,
'L_iterations_per_timestep': True, # save xdmf/h5 data for each LDD iteration for timesteps determined by
# number_of_timesteps_to_analyse. I/O intensive!
'L_iterations_per_timestep': False,
# save solution to xdmf/h5.
'solutions': True, 'solutions': True,
# save absolute differences w.r.t an exact solution to xdmf/h5 file
# to monitor where on the domains errors happen
'absolute_differences': True, 'absolute_differences': True,
# analyise condition numbers for timesteps determined by
# number_of_timesteps_to_analyse and save them over time to csv.
'condition_numbers': analyse_condition, 'condition_numbers': analyse_condition,
# output subsequent iteration errors measured in L^2 to csv for
# timesteps determined by number_of_timesteps_to_analyse.
# Usefull to monitor convergence of the acutal LDD solver.
'subsequent_errors': True 'subsequent_errors': True
} }
else: else:
...@@ -85,9 +119,20 @@ else: ...@@ -85,9 +119,20 @@ else:
'subsequent_errors': True 'subsequent_errors': True
} }
# OUTPUT FILE STRING #########################################################
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
)
##### Domain and Interface #### # DOMAIN AND INTERFACE #######################################################
# global simulation domain domain # global simulation domain domain
sub_domain0_vertices = [df.Point(-1.0, -1.0), sub_domain0_vertices = [df.Point(-1.0, -1.0),
df.Point(1.0, -1.0), df.Point(1.0, -1.0),
...@@ -181,28 +226,19 @@ gravity_acceleration = 9.81 ...@@ -181,28 +226,19 @@ gravity_acceleration = 9.81
L = {# L = {#
# subdom_num : subdomain L for L-scheme # subdom_num : subdomain L for L-scheme
1 : {'wetting' :Lw, 1 : {'wetting' :Lw1,
'nonwetting': Lnw},# 'nonwetting': Lnw1},#
2 : {'wetting' :Lw, 2 : {'wetting' :Lw2,
'nonwetting': Lnw} 'nonwetting': Lnw2}
} }
lambda_param = {# lambda_param = {#
# subdom_num : lambda parameter for the L-scheme # subdom_num : lambda parameter for the L-scheme
1 : {'wetting' :lambda_w, 0 : {'wetting' :lambda_w,
'nonwetting': lambda_nw},# 'nonwetting': lambda_nw},#
2 : {'wetting' :lambda_w,
'nonwetting': lambda_nw}
} }
# intrinsic_permeability = {
# 1: {"wetting": 1,
# "nonwetting": 1},
# 2: {"wetting": 1,
# "nonwetting": 1},
# }
intrinsic_permeability = { intrinsic_permeability = {
1: 1, 1: 1,
2: 1, 2: 1,
...@@ -229,7 +265,7 @@ def rel_perm2w(s): ...@@ -229,7 +265,7 @@ def rel_perm2w(s):
# relative permeabilty wetting on subdomain2 # relative permeabilty wetting on subdomain2
return intrinsic_permeability[2]*s**3 return intrinsic_permeability[2]*s**3
def rel_perm2nw(s): def rel_perm2nw(s):
# relative permeabilty nonwetting on subdosym.cos(0.8*t - (0.8*x + 1/7*y))main2 # relative permeabilty nonwetting on subdomain2
return intrinsic_permeability[2]*(1-s)**3 return intrinsic_permeability[2]*(1-s)**3
_rel_perm2w = ft.partial(rel_perm2w) _rel_perm2w = ft.partial(rel_perm2w)
...@@ -442,7 +478,7 @@ dirichletBC = dict() ...@@ -442,7 +478,7 @@ dirichletBC = dict()
# subdomain index: {outer boudary part index: {phase: expression}} # subdomain index: {outer boudary part index: {phase: expression}}
for subdomain in isRichards.keys(): for subdomain in isRichards.keys():
# if subdomain has no outer boundary, outer_boundary_def_points[subdomain] is None # subdomain can have no outer boundary
if outer_boundary_def_points[subdomain] is None: if outer_boundary_def_points[subdomain] is None:
dirichletBC.update({subdomain: None}) dirichletBC.update({subdomain: None})
else: else:
...@@ -455,13 +491,9 @@ for subdomain in isRichards.keys(): ...@@ -455,13 +491,9 @@ for subdomain in isRichards.keys():
) )
# def saturation(pressure, subdomain_index): # read this file and print it to std out. This way the simulation can produce a
# # inverse capillary pressure-saturation-relationship # log file with ./TP-R-layered_soil.py | tee simulation.log
# return df.conditional(pressure < 0, 1/((1 - pressure)**(1/(subdomain_index + 1))), 1) f = open(thisfile, 'r')
#
# sa
f = open('TP-R-2-patch-mesh-study.py', 'r')
print(f.read()) print(f.read())
f.close() f.close()
...@@ -478,7 +510,8 @@ for starttime in starttimes: ...@@ -478,7 +510,8 @@ for starttime in starttimes:
mesh_study=mesh_study mesh_study=mesh_study
) )
simulation.set_parameters(use_case=use_case, simulation.set_parameters(
use_case=use_case,
output_dir=output_string, output_dir=output_string,
subdomain_def_points=subdomain_def_points, subdomain_def_points=subdomain_def_points,
isRichards=isRichards, isRichards=isRichards,
...@@ -512,9 +545,11 @@ for starttime in starttimes: ...@@ -512,9 +545,11 @@ for starttime in starttimes:
output = simulation.run(analyse_condition=analyse_condition) output = simulation.run(analyse_condition=analyse_condition)
for subdomain_index, subdomain_output in output.items(): for subdomain_index, subdomain_output in output.items():
mesh_h = subdomain_output['mesh_size'] mesh_h = subdomain_output['mesh_size']
for phase, different_errornorms in subdomain_output['errornorm'].items(): for phase, error_dict in subdomain_output['errornorm'].items():
filename = output_dir + "subdomain{}-space-time-errornorm-{}-phase.csv".format(subdomain_index, phase) filename = output_dir \
# for errortype, errornorm in different_errornorms.items(): + "subdomain{}".format(subdomain_index)\
+ "-space-time-errornorm-{}-phase.csv".format(phase)
# for errortype, errornorm in error_dict.items():
# eocfile = open("eoc_filename", "a") # eocfile = open("eoc_filename", "a")
# eocfile.write( str(mesh_h) + " " + str(errornorm) + "\n" ) # eocfile.write( str(mesh_h) + " " + str(errornorm) + "\n" )
...@@ -524,14 +559,25 @@ for starttime in starttimes: ...@@ -524,14 +559,25 @@ for starttime in starttimes:
'mesh_parameter': mesh_resolution, 'mesh_parameter': mesh_resolution,
'mesh_h': mesh_h, 'mesh_h': mesh_h,
} }
for error_type, errornorms in different_errornorms.items(): for norm_type, errornorm in error_dict.items():
data_dict.update( data_dict.update(
{error_type: errornorms} {norm_type: errornorm}
) )
errors = pd.DataFrame(data_dict, index=[mesh_resolution]) errors = pd.DataFrame(data_dict, index=[mesh_resolution])
# check if file exists # check if file exists
if os.path.isfile(filename) == True: if os.path.isfile(filename) is True:
with open(filename, 'a') as f: 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: else:
errors.to_csv(filename, sep='\t', encoding='utf-8', index=False) errors.to_csv(
filename,
sep='\t',
encoding='utf-8',
index=False
)
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment