Select Git revision
TP-R-layered_soil.py
TP-R-layered_soil.py 23.99 KiB
#!/usr/bin/python3
"""Layered 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 functools as ft
# import domainPatch as dp
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")
# init sympy session
sym.init_printing()
# solver_tol = 6E-7
use_case = "TP-R-layered-soil-realistic-new-lambda"
max_iter_num = 600
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
}
# GRID #######################
# mesh_resolution = 20
timestep_size = 0.005
number_of_timesteps = 200
plot_timestep_every = 2
# 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]
Lw1 = 0.25 # /timestep_size
Lnw1 = Lw1
Lw2 = 0.025 # /timestep_size
Lnw2 = Lw2
Lw3 = 0.025 # /timestep_size
Lnw3 = Lw3
Lw4 = 0.025 # /timestep_size
Lnw4 = Lw4
lambda12_w = 4
lambda12_nw = 4
lambda23_w = 40
lambda23_nw = 40
lambda34_w = 40
lambda34_nw = 40
include_gravity = False
debugflag = False
analyse_condition = False
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:
write_to_file = {
'space_errornorms': True,
'meshes_and_markers': True,
'L_iterations_per_timestep': False,
'solutions': True,
'absolute_differences': True,
'condition_numbers': analyse_condition,
'subsequent_errors': True
}
else:
write_to_file = {
'space_errornorms': True,
'meshes_and_markers': True,
'L_iterations_per_timestep': False,
'solutions': True,
'absolute_differences': True,
'condition_numbers': analyse_condition,
'subsequent_errors': True
}
# global domain
subdomain0_vertices = [
df.Point(-1.0, -1.0),
df.Point(1.0, -1.0),
df.Point(1.0, 1.0),
df.Point(-1.0, 1.0)
]
interface12_vertices = [df.Point(-1.0, 0.8),
df.Point(0.3, 0.8),
df.Point(0.5, 0.9),
df.Point(0.8, 0.7),
df.Point(1.0, 0.65)]
# subdomain1.
subdomain1_vertices = [
interface12_vertices[0],
interface12_vertices[1],
interface12_vertices[2],
interface12_vertices[3],
interface12_vertices[4], # southern boundary, 12 interface
subdomain0_vertices[2], # eastern boundary, outer boundary
subdomain0_vertices[3] # northern boundary, outer on_boundary
]
# vertex coordinates of the outer boundaries. If it can not be specified as a
# polygon, use an entry per boundary polygon. This information is used for
# defining the Dirichlet boundary conditions. If a domain is completely
# internal, the dictionary entry should be 0: None
subdomain1_outer_boundary_verts = {
0: [interface12_vertices[4],
subdomain0_vertices[2], # eastern boundary, outer boundary
subdomain0_vertices[3],
interface12_vertices[0]]
}
# interface23
interface23_vertices = [
df.Point(-1.0, 0.0),
df.Point(-0.35, 0.0),
# df.Point(6.5, 4.5),
df.Point(0.0, 0.0),
df.Point(0.5, 0.0),
# df.Point(11.5, 3.5),
# df.Point(13.0, 3)
df.Point(0.85, 0.0),
df.Point(1.0, 0.0)
]
# subdomain1
subdomain2_vertices = [
interface23_vertices[0],
interface23_vertices[1],
interface23_vertices[2],
interface23_vertices[3],
interface23_vertices[4],
interface23_vertices[5], # southern boundary, 23 interface
subdomain1_vertices[4], # eastern boundary, outer boundary
subdomain1_vertices[3],
subdomain1_vertices[2],
subdomain1_vertices[1],
subdomain1_vertices[0] # northern boundary, 12 interface
]
subdomain2_outer_boundary_verts = {
0: [interface23_vertices[5],
subdomain1_vertices[4]],
1: [subdomain1_vertices[0],
interface23_vertices[0]]
}
# interface34
interface34_vertices = [df.Point(-1.0, -0.6),
df.Point(-0.6, -0.45),
df.Point(0.3, -0.25),
df.Point(0.65, -0.6),
df.Point(1.0, -0.7)]
# subdomain3
subdomain3_vertices = [
interface34_vertices[0],
interface34_vertices[1],
interface34_vertices[2],
interface34_vertices[3],
interface34_vertices[4], # southern boundary, 34 interface
subdomain2_vertices[5], # eastern boundary, outer boundary
subdomain2_vertices[4],
subdomain2_vertices[3],
subdomain2_vertices[2],
subdomain2_vertices[1],
subdomain2_vertices[0] # northern boundary, 23 interface
]
subdomain3_outer_boundary_verts = {
0: [interface34_vertices[4],
subdomain2_vertices[5]],
1: [subdomain2_vertices[0],
interface34_vertices[0]]
}
# subdomain4
subdomain4_vertices = [
subdomain0_vertices[0],
subdomain0_vertices[1], # southern boundary, outer boundary
subdomain3_vertices[4], # eastern boundary, outer boundary
subdomain3_vertices[3],
subdomain3_vertices[2],
subdomain3_vertices[1],
subdomain3_vertices[0]
] # northern boundary, 34 interface
subdomain4_outer_boundary_verts = {
0: [subdomain4_vertices[6],
subdomain4_vertices[0],
subdomain4_vertices[1],
subdomain4_vertices[2]]
}
subdomain_def_points = [
subdomain0_vertices,
subdomain1_vertices,
subdomain2_vertices,
subdomain3_vertices,
subdomain4_vertices
]
# interface_vertices introduces a global numbering of interfaces.
interface_def_points = [
interface12_vertices, interface23_vertices, interface34_vertices
]
adjacent_subdomains = [[1, 2], [2, 3], [3, 4]]
# if a subdomain has no outer boundary write None instead, i.e.
# i: None
# if i is the index of the inner subdomain.
outer_boundary_def_points = {
# subdomain number
1: subdomain1_outer_boundary_verts,
2: subdomain2_outer_boundary_verts,
3: subdomain3_outer_boundary_verts,
4: subdomain4_outer_boundary_verts
}
isRichards = {
1: True,
2: True,
3: False,
4: False
}
# isRichards = {
# 1: True,
# 2: True,
# 3: True,
# 4: True
# }
# Dict of the form: { subdom_num : viscosity }
viscosity = {
1: {'wetting': 1,
'nonwetting': 1/50},
2: {'wetting': 1,
'nonwetting': 1/50},
3: {'wetting': 1,
'nonwetting': 1/50},
4: {'wetting': 1,
'nonwetting': 1/50},
}
# 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}}
}
gravity_acceleration = 9.81
# porosities taken from
# https://www.geotechdata.info/parameter/soil-porosity.html
# Dict of the form: { subdom_num : porosity }
porosity = {
1: 0.37, # 0.2, # Clayey gravels, clayey sandy gravels
2: 0.22, # 0.22, # Silty gravels, silty sandy gravels
3: 0.2, # 0.37, # Clayey sands
4: 0.22, # 0.2 # Silty or sandy clay
}
# subdom_num : subdomain L for L-scheme
L = {
1: {'wetting': Lw1,
'nonwetting': Lnw1},
2: {'wetting': Lw2,
'nonwetting': Lnw2},
3: {'wetting': Lw3,
'nonwetting': Lnw3},
4: {'wetting': Lw4,
'nonwetting': Lnw4}
}
# interface_num : lambda parameter for the L-scheme on that interface.
# Note that interfaces are numbered starting from 0, because
# adjacent_subdomains is a list and not a dict. Historic fuckup, I know
lambda_param = {
0: {'wetting': lambda12_w,
'nonwetting': lambda12_nw},
1: {'wetting': lambda23_w,
'nonwetting': lambda23_nw},
2: {'wetting': lambda34_w,
'nonwetting': lambda34_nw},
}
# lambda_param = {
# 1: {'wetting': lambda_w,
# 'nonwetting': lambda_nw},
# 2: {'wetting': lambda_w,
# 'nonwetting': lambda_nw},
# 3: {'wetting': lambda_w,
# 'nonwetting': lambda_nw},
# 4: {'wetting': lambda_w,
# 'nonwetting': lambda_nw},
# }
# after Lewis, see pdf file
intrinsic_permeability = {
1: 1, # sand
2: 0.1, # sand, there is a range
3: 10e-2, # clay has a range
4: 10e-4
}
# relative permeabilty functions on subdomain 1
def rel_perm1w(s):
# relative permeabilty wetting on subdomain1
return intrinsic_permeability[1]*s**2
def rel_perm1nw(s):
# relative permeabilty nonwetting on subdomain1
return intrinsic_permeability[1]*(1-s)**2
# relative permeabilty functions on subdomain 2
def rel_perm2w(s):
# relative permeabilty wetting on subdomain2
return intrinsic_permeability[2]*s**2
def rel_perm2nw(s):
# relative permeabilty nonwetting on subdomain2
return intrinsic_permeability[2]*(1-s)**2
# relative permeabilty functions on subdomain 3
def rel_perm3w(s):
# relative permeabilty wetting on subdomain3
return intrinsic_permeability[3]*s**3
def rel_perm3nw(s):
# relative permeabilty nonwetting on subdomain3
return intrinsic_permeability[3]*(1-s)**3
# relative permeabilty functions on subdomain 4
def rel_perm4w(s):
# relative permeabilty wetting on subdomain4
return intrinsic_permeability[4]*s**3
def rel_perm4nw(s):
# relative permeabilty nonwetting on subdomain4
return intrinsic_permeability[4]*(1-s)**3
_rel_perm1w = ft.partial(rel_perm1w)
_rel_perm1nw = ft.partial(rel_perm1nw)
_rel_perm2w = ft.partial(rel_perm2w)
_rel_perm2nw = ft.partial(rel_perm2nw)
_rel_perm3w = ft.partial(rel_perm3w)
_rel_perm3nw = ft.partial(rel_perm3nw)
_rel_perm4w = ft.partial(rel_perm4w)
_rel_perm4nw = ft.partial(rel_perm4nw)
subdomain1_rel_perm = {
'wetting': _rel_perm1w,
'nonwetting': _rel_perm1nw
}
subdomain2_rel_perm = {
'wetting': _rel_perm2w,
'nonwetting': _rel_perm2nw
}
subdomain3_rel_perm = {
'wetting': _rel_perm3w,
'nonwetting': _rel_perm3nw
}
subdomain4_rel_perm = {
'wetting': _rel_perm4w,
'nonwetting': _rel_perm4nw
}
# dictionary of relative permeabilties on all domains.
# relative_permeability = {
# 1: subdomain1_rel_perm,
# 2: subdomain1_rel_perm,
# 3: subdomain2_rel_perm,
# 4: subdomain2_rel_perm
# }
relative_permeability = {
1: subdomain1_rel_perm,
2: subdomain2_rel_perm,
3: subdomain3_rel_perm,
4: subdomain4_rel_perm
}
# definition of the derivatives of the relative permeabilities
# relative permeabilty functions on subdomain 1
def rel_perm1w_prime(s):
# relative permeabilty on subdomain1
return intrinsic_permeability[1]*2*s
def rel_perm1nw_prime(s):
# relative permeabilty on subdomain1
return -1*intrinsic_permeability[1]*2*(1-s)
def rel_perm2w_prime(s):
# relative permeabilty on subdomain2
return intrinsic_permeability[2]*2*s
def rel_perm2nw_prime(s):
# relative permeabilty on subdomain2
return -1*intrinsic_permeability[2]*2*(1-s)
# definition of the derivatives of the relative permeabilities
# relative permeabilty functions on subdomain 3
def rel_perm3w_prime(s):
# relative permeabilty on subdomain3
return intrinsic_permeability[3]*3*s**2
def rel_perm3nw_prime(s):
# relative permeabilty on subdomain3
return -1*intrinsic_permeability[3]*3*(1-s)**2
# definition of the derivatives of the relative permeabilities
# relative permeabilty functions on subdomain 4
def rel_perm4w_prime(s):
# relative permeabilty on subdomain4
return intrinsic_permeability[4]*3*s**2
def rel_perm4nw_prime(s):
# relative permeabilty on subdomain4
return -1*intrinsic_permeability[4]*3*(1-s)**2
_rel_perm1w_prime = ft.partial(rel_perm1w_prime)
_rel_perm1nw_prime = ft.partial(rel_perm1nw_prime)
_rel_perm2w_prime = ft.partial(rel_perm2w_prime)
_rel_perm2nw_prime = ft.partial(rel_perm2nw_prime)
_rel_perm3w_prime = ft.partial(rel_perm3w_prime)
_rel_perm3nw_prime = ft.partial(rel_perm3nw_prime)
_rel_perm4w_prime = ft.partial(rel_perm4w_prime)
_rel_perm4nw_prime = ft.partial(rel_perm4nw_prime)
subdomain1_rel_perm_prime = {
'wetting': _rel_perm1w_prime,
'nonwetting': _rel_perm1nw_prime
}
subdomain2_rel_perm_prime = {
'wetting': _rel_perm2w_prime,
'nonwetting': _rel_perm2nw_prime
}
subdomain3_rel_perm_prime = {
'wetting': _rel_perm3w_prime,
'nonwetting': _rel_perm3nw_prime
}
subdomain4_rel_perm_prime = {
'wetting': _rel_perm4w_prime,
'nonwetting': _rel_perm4nw_prime
}
# dictionary of relative permeabilties on all domains.
# ka_prime = {
# 1: subdomain1_rel_perm_prime,
# 2: subdomain1_rel_perm_prime,
# 3: subdomain2_rel_perm_prime,
# 4: subdomain2_rel_perm_prime
# }
ka_prime = {
1: subdomain1_rel_perm_prime,
2: subdomain2_rel_perm_prime,
3: subdomain3_rel_perm_prime,
4: subdomain4_rel_perm_prime
}
# 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
# this function needs to be monotonically decreasing in the capillary pressure
# pc. Since in the richards case pc=-pw, this becomes as a function of pw a
# mono tonically 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(pc, n_index, alpha):
# inverse capillary pressure-saturation-relationship
expr = df.conditional(
pc > 0, 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index)), 1)
return expr
# 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
def saturation_sym(pc, n_index, alpha):
# inverse capillary pressure-saturation-relationship
# df.conditional(pc > 0,
return 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index))
# derivative of S-pc relationship with respect to pc. This is needed for the
# construction of a analytic solution.
def saturation_sym_prime(pc, n_index, alpha):
# inverse capillary pressure-saturation-relationship
expr = -(alpha*(n_index - 1)*(alpha*pc)**(n_index - 1))\
/ ((1 + (alpha*pc)**n_index)**((2*n_index - 1)/n_index))
return expr
# note that the conditional definition of S-pc in the nonsymbolic part will be
# incorporated in the construction of the exact solution below.
S_pc_sym = {
1: ft.partial(saturation_sym, n_index=3, alpha=0.001),
2: ft.partial(saturation_sym, n_index=3, alpha=0.001),
3: ft.partial(saturation_sym, n_index=6, alpha=0.001),
4: ft.partial(saturation_sym, n_index=6, alpha=0.001)
}
S_pc_sym_prime = {
1: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
2: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
3: ft.partial(saturation_sym_prime, n_index=6, alpha=0.001),
4: ft.partial(saturation_sym_prime, n_index=6, alpha=0.001)
}
sat_pressure_relationship = {
1: ft.partial(saturation, n_index=3, alpha=0.001),
2: ft.partial(saturation, n_index=3, alpha=0.001),
3: ft.partial(saturation, n_index=6, alpha=0.001),
4: ft.partial(saturation, n_index=6, alpha=0.001)
}
#############################################
# Manufacture source expressions with sympy #
#############################################
x, y = sym.symbols('x[0], x[1]') # needed by UFL
t = sym.symbols('t', positive=True)
p_e_sym_2patch = {
1: {'wetting': -6 - (1+t*t)*(1 + x*x + y*y),
'nonwetting': 0.0*t}, # -1-t*(1.1 + y + x**2)**2},
2: {'wetting': -6.0 - (1.0 + t*t)*(1.0 + x*x),
'nonwetting': (-1-t*(1.1+y + x**2))*y**2},
# 'nonwetting': (-1-t*(1.1 + x**2)**2 - sym.sqrt(5+t**2))*y**2},
}
p_e_sym = {
1: {'wetting': p_e_sym_2patch[1]['wetting'],
'nonwetting': p_e_sym_2patch[1]['nonwetting']},
2: {'wetting': p_e_sym_2patch[1]['wetting'],
'nonwetting': p_e_sym_2patch[1]['nonwetting']},
3: {'wetting': p_e_sym_2patch[2]['wetting'],
'nonwetting': p_e_sym_2patch[2]['nonwetting']},
4: {'wetting': p_e_sym_2patch[2]['wetting'],
'nonwetting': p_e_sym_2patch[2]['nonwetting']}
}
# p_e_sym = {
# 1: {'wetting': 1.0 - (1.0 + t*t)*(10.0 + x*x + (y-5.0)*(y-5.0)),
# 'nonwetting':
# -2 - t*(1 + (y - 5.0) + x**2)**2
# - sym.sqrt(2 + t**2)*(1 + (y - 5.0))
# },
# 2: {'wetting': 1.0 - (1.0 + t*t)*(10.0 + x*x + (y-5.0)*(y-5.0)),
# 'nonwetting':
# - 2 - t*(1 + (y - 5.0) + x**2)**2
# - sym.sqrt(2 + t**2)*(1 + (y - 5.0))
# },
# 3: {'wetting':
# 1.0 - (1.0 + t*t)*(10.0 + x*x + (y - 5.0)*(y - 5.0))
# - (y - 5.0)*(y - 5.0)*3*sym.sin(-2*t + 2*x)*sym.sin(1/2*y - 1.2*t),
# 'nonwetting': - 2 - t*(1 + x**2)**2 - sym.sqrt(2 + t**2)
# },
# 4: {'wetting':
# 1.0 - (1.0 + t*t)*(10.0 + x*x + (y - 5.0)*(y - 5.0))
# - (y - 5.0)*(y - 5.0)*3*sym.sin(-2*t + 2*x)*sym.sin(1/2*y - 1.2*t),
# 'nonwetting': - 2 - t*(1 + x**2)**2 - sym.sqrt(2 + t**2)
# }
# }
pc_e_sym = dict()
for subdomain, isR in isRichards.items():
if isR:
pc_e_sym.update({subdomain: -p_e_sym[subdomain]['wetting']})
else:
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,
)
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()
# 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}}
for subdomain in isRichards.keys():
# if subdomain has no outer boundary, outer_boundary_def_points[subdomain]
# is None
if outer_boundary_def_points[subdomain] is None:
dirichletBC.update({subdomain: None})
else:
dirichletBC.update({subdomain: dict()})
# set the dirichlet conditions to be the same code as exact solution on
# the subdomain.
for outer_boundary_ind in outer_boundary_def_points[subdomain].keys():
dirichletBC[subdomain].update(
{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')
print(f.read())
f.close()
for starttime in starttimes:
for mesh_resolution, solver_tol in resolutions.items():
# initialise LDD simulation class
simulation = ldd.LDDsimulation(
tol=1E-14,
LDDsolver_tol=solver_tol,
debug=debugflag,
max_iter_num=max_iter_num,
FEM_Lagrange_degree=FEM_Lagrange_degree,
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.initialise()
output_dir = simulation.output_dir
# simulation.write_exact_solution_to_xdmf()
output = simulation.run(analyse_condition=analyse_condition)
for subdomain_index, subdomain_output in output.items():
mesh_h = subdomain_output['mesh_size']
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 norm_type, errornorm in error_dict.items():
data_dict.update(
{norm_type: errornorm}
)
errors = pd.DataFrame(data_dict, index=[mesh_resolution])
# check if file exists
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
)
else:
errors.to_csv(
filename,
sep='\t',
encoding='utf-8',
index=False
)