From 84518db196fe11febaefa0aa6f0485db9a4faa8c Mon Sep 17 00:00:00 2001
From: David <forenkram@gmx.de>
Date: Tue, 9 Jun 2020 11:13:27 +0200
Subject: [PATCH] update and clean up example
---
.../mesh_studies/TP-R-2-patch-mesh-study.py | 228 +++++++++++-------
1 file changed, 137 insertions(+), 91 deletions(-)
diff --git a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study.py
index 1feb5db..6afe2dc 100755
--- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study.py
+++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study.py
@@ -1,28 +1,45 @@
#!/usr/bin/python3
+"""TPR 2 patch soil simulation.
+
+This program sets up an LDD simulation
+"""
+
import dolfin as df
-import mshr
-import numpy as np
import sympy as sym
-import typing as tp
-import domainPatch as dp
-import LDDsimulation as ldd
import functools as ft
+import LDDsimulation as ldd
import helpers as hlp
import datetime
import os
import pandas as pd
-date = datetime.datetime.now()
-datestr = date.strftime("%Y-%m-%d")
-#import ufl as ufl
-
# init sympy session
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"
-# 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
FEM_Lagrange_degree = 1
+
+# GRID AND MESH STUDY SPECIFICATIONS #########################################
mesh_study = True
resolutions = {
1: 1e-6,
@@ -36,18 +53,19 @@ resolutions = {
256: 1e-6,
}
-############ GRID #######################
-# mesh_resolution = 20
+# starttimes gives a list of starttimes to run the simulation from.
+# 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
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
-Lnw= 0.025
+# LDD scheme parameters ######################################################
+Lw1 = 0.025 #/timestep_size
+Lnw1= 0.025
+
+Lw2 = 0.025 #/timestep_size
+Lnw2= 0.025
lambda_w = 40
lambda_nw = 40
@@ -56,22 +74,38 @@ include_gravity = False
debugflag = False
analyse_condition = True
-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)
+# I/O CONFIG #################################################################
+# when number_of_timesteps is high, it might take a long time to write all
+# timesteps to disk. Therefore, you can choose to only write data of every
+# plot_timestep_every timestep to disk.
+plot_timestep_every = 4
+# 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:
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,
+ # save the mesh and marker functions to disk
'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,
+ # save absolute differences w.r.t an exact solution to xdmf/h5 file
+ # to monitor where on the domains errors happen
'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,
+ # 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
}
else:
@@ -85,9 +119,20 @@ else:
'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
sub_domain0_vertices = [df.Point(-1.0, -1.0),
df.Point(1.0, -1.0),
@@ -181,28 +226,19 @@ gravity_acceleration = 9.81
L = {#
# subdom_num : subdomain L for L-scheme
- 1 : {'wetting' :Lw,
- 'nonwetting': Lnw},#
- 2 : {'wetting' :Lw,
- 'nonwetting': Lnw}
+ 1 : {'wetting' :Lw1,
+ 'nonwetting': Lnw1},#
+ 2 : {'wetting' :Lw2,
+ 'nonwetting': Lnw2}
}
lambda_param = {#
# subdom_num : lambda parameter for the L-scheme
- 1 : {'wetting' :lambda_w,
+ 0 : {'wetting' :lambda_w,
'nonwetting': lambda_nw},#
- 2 : {'wetting' :lambda_w,
- 'nonwetting': lambda_nw}
}
-# intrinsic_permeability = {
-# 1: {"wetting": 1,
-# "nonwetting": 1},
-# 2: {"wetting": 1,
-# "nonwetting": 1},
-# }
-
intrinsic_permeability = {
1: 1,
2: 1,
@@ -229,7 +265,7 @@ def rel_perm2w(s):
# relative permeabilty wetting on subdomain2
return intrinsic_permeability[2]*s**3
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
_rel_perm2w = ft.partial(rel_perm2w)
@@ -442,7 +478,7 @@ dirichletBC = dict()
# subdomain index: {outer boudary part index: {phase: expression}}
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:
dirichletBC.update({subdomain: None})
else:
@@ -455,13 +491,9 @@ for subdomain in isRichards.keys():
)
-# def saturation(pressure, subdomain_index):
-# # inverse capillary pressure-saturation-relationship
-# return df.conditional(pressure < 0, 1/((1 - pressure)**(1/(subdomain_index + 1))), 1)
-#
-# sa
-
-f = open('TP-R-2-patch-mesh-study.py', 'r')
+# 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(thisfile, 'r')
print(f.read())
f.close()
@@ -478,33 +510,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
@@ -512,26 +545,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
+ )
--
GitLab