From 9d289a6aed3b5e690555293c5e160c612afe4f3e Mon Sep 17 00:00:00 2001
From: David Seus <david.seus@ians.uni-stuttgart.de>
Date: Tue, 24 Sep 2019 15:21:36 +0200
Subject: [PATCH] set up TP-TP layered soil with inner patch
---
...nner_patch-realistic-split-up-interface.py | 750 ------------------
...layered_soil_with_inner_patch-realistic.py | 205 +++--
2 files changed, 145 insertions(+), 810 deletions(-)
delete mode 100755 Two-phase-Two-phase/multi-patch/TP-TP-layered-soil-case-with-inner-patch/TP-TP-layered_soil_with_inner_patch-realistic-split-up-interface.py
diff --git a/Two-phase-Two-phase/multi-patch/TP-TP-layered-soil-case-with-inner-patch/TP-TP-layered_soil_with_inner_patch-realistic-split-up-interface.py b/Two-phase-Two-phase/multi-patch/TP-TP-layered-soil-case-with-inner-patch/TP-TP-layered_soil_with_inner_patch-realistic-split-up-interface.py
deleted file mode 100755
index bf29769..0000000
--- a/Two-phase-Two-phase/multi-patch/TP-TP-layered-soil-case-with-inner-patch/TP-TP-layered_soil_with_inner_patch-realistic-split-up-interface.py
+++ /dev/null
@@ -1,750 +0,0 @@
-#!/usr/bin/python3
-"""This program sets up a domain together with a decomposition into subdomains
-modelling layered soil. This is used for our LDD article with tp-tp and tp-r
-coupling.
-
-Along with the subdomains and the mesh domain markers are set upself.
-The resulting mesh is saved into files for later use.
-"""
-
-#!/usr/bin/python3
-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
-
-# init sympy session
-sym.init_printing()
-
-use_case = "TP-TP-layered-soil-with-inner-patch-realistic-split-interface-with-gravity"
-solver_tol = 5E-6
-max_iter_num = 100
-############ GRID #######################ü
-mesh_resolution = 30
-timestep_size = 0.00001
-number_of_timesteps = 1
-# 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 = 0
-starttime = 0
-
-Lw = 2 #/timestep_size
-Lnw=Lw
-
-lambda_w = 40
-lambda_nw = 40
-
-include_gravity = True
-debugflag = True
-analyse_condition = False
-
-output_string = "./output/2019-08-30-debug-{}-timesteps{}_".format(use_case, number_of_timesteps)
-
-# 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)]
-
-
- # 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)]
-
-interface24_vertices = [interface23_vertices[2],
- df.Point(0.6, 0.0),
- ]
-
-interface25_vertices = [interface24_vertices[1],
- df.Point(1.0, 0.0)
- ]
-
-
-interface32_vertices = [interface23_vertices[2],
- interface23_vertices[1],
- interface23_vertices[0]]
-
-
-interface36_vertices = [df.Point(-1.0, -0.6),
- df.Point(-0.6, -0.45)]
-
-
-interface46_vertices = [interface36_vertices[1],
- df.Point(0.3, -0.25)]
-
-interface56_vertices = [interface46_vertices[1],
- df.Point(0.65, -0.6),
- df.Point(1.0, -0.7)]
-
-
-
-
-interface34_vertices = [interface36_vertices[1],
- interface23_vertices[2]]
-# interface36
-
-
-interface45_vertices_a = [interface56_vertices[0],
- df.Point(0.7, -0.2),#df.Point(0.7, -0.2),
- ]
-interface45_vertices_b = [df.Point(0.7, -0.2),#df.Point(0.7, -0.2),
- interface25_vertices[0]
- ]
-
-
-# # 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]]
-# }
-#
-
-
-# #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]]
-# }
-#
-
-# interface_vertices introduces a global numbering of interfaces.
-interface_def_points = [interface12_vertices,
- interface23_vertices,
- interface24_vertices,
- interface25_vertices,
- interface34_vertices,
- interface36_vertices,
- # interface45_vertices,
- interface45_vertices_a,
- interface45_vertices_b,
- interface46_vertices,
- interface56_vertices,
- ]
-adjacent_subdomains = [[1,2],
- [2,3],
- [2,4],
- [2,5],
- [3,4],
- [3,6],
- [4,5],
- [4,5],
- [4,6],
- [5,6]
- ]
-
-# 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: [subdomain1_vertices[4], #
- subdomain1_vertices[5], # eastern boundary, outer boundary
- subdomain1_vertices[6],
- subdomain1_vertices[0]]
-}
-
-#subdomain1
-subdomain2_vertices = [interface23_vertices[0],
- interface23_vertices[1],
- interface23_vertices[2],
- interface24_vertices[1],
- interface25_vertices[1], # 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: [subdomain2_vertices[9],
- subdomain2_vertices[0]],
- 1: [subdomain2_vertices[4],
- subdomain2_vertices[5]]
-}
-
-
-subdomain3_vertices = [interface36_vertices[0],
- interface36_vertices[1],
- # interface34_vertices[0],
- interface34_vertices[1],
- # interface32_vertices[0],
- interface32_vertices[1],
- interface32_vertices[2]
- ]
-
-subdomain3_outer_boundary_verts = {
- 0: [subdomain3_vertices[4],
- subdomain3_vertices[0]]
-}
-
-
-# subdomain3
-subdomain4_vertices = [interface46_vertices[0],
- interface46_vertices[1],
- # interface45_vertices[1],
- interface45_vertices_a[1],
- interface24_vertices[1],
- interface24_vertices[0],
- interface34_vertices[1]
- ]
-
-subdomain4_outer_boundary_verts = None
-
-subdomain5_vertices = [interface56_vertices[0],
- interface56_vertices[1],
- interface56_vertices[2],
- interface25_vertices[1],
- interface25_vertices[0],
- interface45_vertices_b[1],
- interface45_vertices_b[0]
-]
-
-
-subdomain5_outer_boundary_verts = {
- 0: [subdomain5_vertices[2],
- subdomain5_vertices[3]]
-}
-
-
-
-subdomain6_vertices = [subdomain0_vertices[0],
- subdomain0_vertices[1], # southern boundary, outer boundary
- interface56_vertices[2],
- interface56_vertices[1],
- interface56_vertices[0],
- interface36_vertices[1],
- interface36_vertices[0]
- ]
-
-subdomain6_outer_boundary_verts = {
- 0: [subdomain6_vertices[6],
- subdomain6_vertices[0],
- subdomain6_vertices[1],
- subdomain6_vertices[2]]
-}
-
-
-subdomain_def_points = [subdomain0_vertices,#
- subdomain1_vertices,#
- subdomain2_vertices,#
- subdomain3_vertices,#
- subdomain4_vertices,
- subdomain5_vertices,
- subdomain6_vertices
- ]
-
-
-# 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,
- 5: subdomain5_outer_boundary_verts,
- 6: subdomain6_outer_boundary_verts
-}
-
-
-isRichards = {
- 1: False,
- 2: False,
- 3: False,
- 4: False,
- 5: False,
- 6: False
- }
-
-# isRichards = {
-# 1: True,
-# 2: True,
-# 3: True,
-# 4: True,
-# 5: True,
-# 6: 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},
- 5: {'wetting' :1,
- 'nonwetting': 1/50},
- 6: {'wetting' :1,
- 'nonwetting': 1/50},
-}
-
-# Dict of the form: { subdom_num : density }
-densities = {
- 1: {'wetting': 997, #997
- 'nonwetting': 1.225}, #1}, #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}
- 5: {'wetting': 997, #997
- 'nonwetting': 1.225}, #1.225},
- 6: {'wetting': 997, #997
- 'nonwetting': 1.225} #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.2, #0.2, # Clayey gravels, clayey sandy gravels
- 2: 0.22, #0.22, # Silty gravels, silty sandy gravels
- 3: 0.22, #0.37, # Clayey sands
- 4: 0.27, #0.2 # Silty or sandy clay
- 5: 0.2, #
- 6: 0.02, #
-}
-
-# subdom_num : subdomain L for L-scheme
-L = {
- 1: {'wetting' :Lw,
- 'nonwetting': Lnw},
- 2: {'wetting' :Lw,
- 'nonwetting': Lnw},
- 3: {'wetting' :Lw,
- 'nonwetting': Lnw},
- 4: {'wetting' :Lw,
- 'nonwetting': Lnw},
- 5: {'wetting' :Lw,
- 'nonwetting': Lnw},
- 6: {'wetting' :Lw,
- 'nonwetting': Lnw}
-}
-
-# subdom_num : lambda parameter for the L-scheme
-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},#
- 5: {'wetting': lambda_w,
- 'nonwetting': lambda_nw},#
- 6: {'wetting': lambda_w,
- 'nonwetting': lambda_nw},#
-}
-
-
-## relative permeabilty functions on subdomain 1
-def rel_perm1w(s):
- # relative permeabilty wetting on subdomain1
- return s**2
-
-
-def rel_perm1nw(s):
- # relative permeabilty nonwetting on subdomain1
- return (1-s)**2
-
-
-## relative permeabilty functions on subdomain 2
-def rel_perm2w(s):
- # relative permeabilty wetting on subdomain2
- return s**3
-
-
-def rel_perm2nw(s):
- # relative permeabilty nonwetting on subdosym.cos(0.8*t - (0.8*x + 1/7*y))main2
- return (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)
-
-subdomain1_rel_perm = {
- 'wetting': _rel_perm1w,#
- 'nonwetting': _rel_perm1nw
-}
-
-subdomain2_rel_perm = {
- 'wetting': _rel_perm2w,#
- 'nonwetting': _rel_perm2nw
-}
-
-# _rel_perm3 = ft.partial(rel_perm2)
-# subdomain3_rel_perm = subdomain2_rel_perm.copy()
-#
-# _rel_perm4 = ft.partial(rel_perm1)
-# subdomain4_rel_perm = subdomain1_rel_perm.copy()
-
-# 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,
- 5: subdomain2_rel_perm,
- 6: subdomain2_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 2*s
-
-def rel_perm1nw_prime(s):
- # relative permeabilty on subdomain1
- return -2*(1-s)
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 1
-def rel_perm2w_prime(s):
- # relative permeabilty on subdomain1
- return 3*s**2
-
-def rel_perm2nw_prime(s):
- # relative permeabilty on subdomain1
- return -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)
-
-subdomain1_rel_perm_prime = {
- 'wetting': _rel_perm1w_prime,
- 'nonwetting': _rel_perm1nw_prime
-}
-
-
-subdomain2_rel_perm_prime = {
- 'wetting': _rel_perm2w_prime,
- 'nonwetting': _rel_perm2nw_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,
- 5: subdomain2_rel_perm_prime,
- 6: subdomain2_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
-# return df.conditional(pc > 0, 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index)), 1)
-#
-# # 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
-# return -(alpha*(n_index - 1)*(alpha*pc)**(n_index - 1)) / ( (1 + (alpha*pc)**n_index)**((2*n_index - 1)/n_index) )
-#
-# derivative of S-pc relationship with respect to pc. This is needed for the
-# construction of a analytic solution.
-
-#
-# # 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=3, alpha=0.001),
-# 4: ft.partial(saturation_sym, n_index=3, alpha=0.001),
-# 5: ft.partial(saturation_sym, n_index=3, alpha=0.001),
-# 6: ft.partial(saturation_sym, n_index=3, 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=3, alpha=0.001),
-# 4: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
-# 5: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
-# 6: ft.partial(saturation_sym_prime, n_index=3, 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=3, alpha=0.001),
-# 4: ft.partial(saturation, n_index=3, alpha=0.001),
-# 5: ft.partial(saturation, n_index=3, alpha=0.001),
-# 6: ft.partial(saturation, n_index=3, alpha=0.001)
-# }
-
-def saturation(pc, n_index):
- # inverse capillary pressure-saturation-relationship
- return df.conditional(pc > 0, 1/((1 + pc)**(1/(n_index + 1))), 1)
-
-
-def saturation_sym(pc, n_index):
- # inverse capillary pressure-saturation-relationship
- return 1/((1 + pc)**(1/(n_index + 1)))
-
-def saturation_sym_prime(pc, n_index):
- # inverse capillary pressure-saturation-relationship
- return -1/((n_index+1)*(1 + pc)**((n_index+2)/(n_index+1)))
-
-
-S_pc_sym = {
- 1: ft.partial(saturation_sym, n_index=1),
- 2: ft.partial(saturation_sym, n_index=1),
- 3: ft.partial(saturation_sym, n_index=2),
- 4: ft.partial(saturation_sym, n_index=2),
- 5: ft.partial(saturation_sym, n_index=2),
- 6: ft.partial(saturation_sym, n_index=2)
-}
-
-S_pc_sym_prime = {
- 1: ft.partial(saturation_sym_prime, n_index=1),
- 2: ft.partial(saturation_sym_prime, n_index=1),
- 3: ft.partial(saturation_sym_prime, n_index=2),
- 4: ft.partial(saturation_sym_prime, n_index=2),
- 5: ft.partial(saturation_sym_prime, n_index=2),
- 6: ft.partial(saturation_sym_prime, n_index=2)
-}
-
-sat_pressure_relationship = {
- 1: ft.partial(saturation, n_index=1),
- 2: ft.partial(saturation, n_index=1),
- 3: ft.partial(saturation, n_index=2),
- 4: ft.partial(saturation, n_index=2),
- 5: ft.partial(saturation, n_index=2),
- 6: ft.partial(saturation, n_index=2)
-}
-
-
-#############################################
-# 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 = {
- 1: {'wetting': -5.0 - (1.0 + t*t)*(1.0 + x*x + y*y),
- 'nonwetting': (-1 -t*(1.1 + y + x**2)) },
- 2: {'wetting': -5.0 - (1.0 + t*t)*(1.0 + x*x + y*y),
- 'nonwetting': (-1 -t*(1.1 + y + x**2)) },
- 3: {'wetting': (-5.0 - (1.0 + t*t)*(1.0 + x*x)),
- 'nonwetting': (-1 -t*(1 + x**2) - sym.sqrt(2+t**2)*(1+y)*y**2) },
- 4: {'wetting': (-5.0 - (1.0 + t*t)*(1.0 + x*x)),
- 'nonwetting': (-1 -t*(1 + x**2) - sym.sqrt(2+t**2)*(1+y)*y**2) },
- 5: {'wetting': (-5.0 - (1.0 + t*t)*(1.0 + x*x)),
- 'nonwetting': (-1 -t*(1 + x**2) - sym.sqrt(2+t**2)*(1+y)*y**2) },
- 6: {'wetting': (-5.0 - (1.0 + t*t)*(1.0 + x*x)),
- 'nonwetting': (-1 -t*(1 + x**2) - sym.sqrt(2+t**2)*(1+y)*y**2) },
- # 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)*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)*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 = {
-# 1: p_e_sym[1]['nonwetting'] - p_e_sym[1]['wetting'],
-# 2: p_e_sym[2]['nonwetting'] - p_e_sym[2]['wetting'],
-# 3: p_e_sym[3]['nonwetting'] - p_e_sym[3]['wetting'],
-# 4: p_e_sym[4]['nonwetting'] - p_e_sym[4]['wetting'],
-# 5: p_e_sym[5]['nonwetting'] - p_e_sym[5]['wetting'],
-# 6: p_e_sym[5]['nonwetting'] - p_e_sym[6]['wetting']
-# }
-
-# pc_e_sym = {
-# 1: -p_e_sym[1]['wetting'],
-# 2: -p_e_sym[2]['wetting'],
-# 3: -p_e_sym[3]['wetting'],
-# 4: -p_e_sym[4]['wetting'],
-# 5: -p_e_sym[5]['wetting'],
-# 6: -p_e_sym[6]['wetting']
-# }
-
-
-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]}
- )
-
-write_to_file = {
- 'meshes_and_markers': True,
- 'L_iterations': True
-}
-
-# initialise LDD simulation class
-simulation = ldd.LDDsimulation(
- tol=1E-14,
- debug=debugflag,
- LDDsolver_tol=solver_tol,
- max_iter_num=max_iter_num
- )
-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,
- 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()
-# print(simulation.__dict__)
-simulation.run(analyse_condition=analyse_condition)
-# simulation.LDDsolver(time=0, debug=True, analyse_timestep=True)
-# df.info(parameters, True)
diff --git a/Two-phase-Two-phase/multi-patch/TP-TP-layered-soil-case-with-inner-patch/TP-TP-layered_soil_with_inner_patch-realistic.py b/Two-phase-Two-phase/multi-patch/TP-TP-layered-soil-case-with-inner-patch/TP-TP-layered_soil_with_inner_patch-realistic.py
index 37d337a..ec01e72 100755
--- a/Two-phase-Two-phase/multi-patch/TP-TP-layered-soil-case-with-inner-patch/TP-TP-layered_soil_with_inner_patch-realistic.py
+++ b/Two-phase-Two-phase/multi-patch/TP-TP-layered-soil-case-with-inner-patch/TP-TP-layered_soil_with_inner_patch-realistic.py
@@ -17,33 +17,79 @@ 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()
-
-use_case = "TP-TP-layered-soil-with-inner-patch-realistic-with-gravity"
-solver_tol = 2E-6
-
-############ GRID #######################ü
-mesh_resolution = 20
-timestep_size = 0.00005
-number_of_timesteps = 1000
+# solver_tol = 6E-7
+use_case = "TP-TP-layered-soil-realistic"
+max_iter_num = 50
+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: 5e-7, # h=0.1412
+ # 32: 4e-7, # h=0.0706
+ 64: 7e-7,
+ # 128: 5e-7
+ }
+
+############ GRID #######################
+# mesh_resolution = 20
+timestep_size = 0.005
+number_of_timesteps = 15
+plot_timestep_every = 1
# 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
+number_of_timesteps_to_analyse = 0
+starttime = 0.0
-Lw = 1 #/timestep_size
+Lw = 0.025 #/timestep_size
Lnw=Lw
-lambda_w = 4
-lambda_nw = 4
+lambda_w = 40
+lambda_nw = 40
-include_gravity = True
+include_gravity = False
debugflag = True
analyse_condition = False
-output_string = "./output/2019-08-30-{}_timesteps{}_".format(use_case, number_of_timesteps)
+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': False,
+ 'absolute_differences': False,
+ 'condition_numbers': analyse_condition,
+ 'subsequent_errors': False
+ }
+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), #
@@ -96,11 +142,15 @@ interface34_vertices = [interface36_vertices[1],
interface23_vertices[2]]
# interface36
-interface45_vertices = [interface56_vertices[0],
- df.Point(0.7, -0.2),
+
+interface45_vertices_a = [interface56_vertices[0],
+ df.Point(0.7, -0.2),#df.Point(0.7, -0.2),
+ ]
+interface45_vertices_b = [df.Point(0.7, -0.2),#df.Point(0.7, -0.2),
interface25_vertices[0]
]
+
# # subdomain1.
# subdomain1_vertices = [interface12_vertices[0],
# interface12_vertices[1],
@@ -151,7 +201,9 @@ interface_def_points = [interface12_vertices,
interface25_vertices,
interface34_vertices,
interface36_vertices,
- interface45_vertices,
+ # interface45_vertices,
+ interface45_vertices_a,
+ interface45_vertices_b,
interface46_vertices,
interface56_vertices,
]
@@ -162,6 +214,7 @@ adjacent_subdomains = [[1,2],
[3,4],
[3,6],
[4,5],
+ [4,5],
[4,6],
[5,6]
]
@@ -224,7 +277,8 @@ subdomain3_outer_boundary_verts = {
# subdomain3
subdomain4_vertices = [interface46_vertices[0],
interface46_vertices[1],
- interface45_vertices[1],
+ # interface45_vertices[1],
+ interface45_vertices_a[1],
interface24_vertices[1],
interface24_vertices[0],
interface34_vertices[1]
@@ -237,10 +291,11 @@ subdomain5_vertices = [interface56_vertices[0],
interface56_vertices[2],
interface25_vertices[1],
interface25_vertices[0],
- interface45_vertices[1],
- interface45_vertices[0]
+ interface45_vertices_b[1],
+ interface45_vertices_b[0]
]
+
subdomain5_outer_boundary_verts = {
0: [subdomain5_vertices[2],
subdomain5_vertices[3]]
@@ -345,7 +400,7 @@ gravity_acceleration = 9.81
# Dict of the form: { subdom_num : porosity }
porosity = {
1: 0.2, #0.2, # Clayey gravels, clayey sandy gravels
- 2: 0.22, #0.22, # Silty gravels, silty sandy gravels
+ 2: 0.0022, #0.22, # Silty gravels, silty sandy gravels
3: 0.22, #0.37, # Clayey sands
4: 0.27, #0.2 # Silty or sandy clay
5: 0.2, #
@@ -695,42 +750,72 @@ for subdomain in isRichards.keys():
{outer_boundary_ind: exact_solution[subdomain]}
)
-write_to_file = {
- 'meshes_and_markers': True,
- 'L_iterations': True
-}
-# initialise LDD simulation class
-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,
- 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,
- 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()
-# print(simulation.__dict__)
-simulation.run(analyse_condition=analyse_condition)
-# simulation.LDDsolver(time=0, debug=True, analyse_timestep=True)
-# df.info(parameters, True)
+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, 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
+ data_dict = {
+ 'mesh_parameter': mesh_resolution,
+ 'mesh_h': mesh_h,
+ }
+ for error_type, errornorms in different_errornorms.items():
+ data_dict.update(
+ {error_type: errornorms}
+ )
+ errors = pd.DataFrame(data_dict, index=[mesh_resolution])
+ # check if file exists
+ if os.path.isfile(filename) == 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)
--
GitLab