From a6c6b8bf3e752d5e5d0f2655dead3d90324b1e71 Mon Sep 17 00:00:00 2001
From: David Seus <david.seus@ians.uni-stuttgart.de>
Date: Tue, 13 Oct 2020 13:30:40 +0200
Subject: [PATCH] adjust readme, bring TP-R-2-patch-test.py up to par.

---
 Usecases/README.md                            | 255 +++++++++++++++++-
 .../TP-R-2-patch-test.py                      | 183 ++-----------
 2 files changed, 271 insertions(+), 167 deletions(-)

diff --git a/Usecases/README.md b/Usecases/README.md
index d3d9b0e..5fcb5c9 100644
--- a/Usecases/README.md
+++ b/Usecases/README.md
@@ -68,4 +68,257 @@ simulation data. Chose a descriptive name for the usecase you are simulating.!
 - `thisfile` should be the name of the script file. This variable is used to
 output the content of the script to stdout and into the logfile if the script
 `run_simulation` is used to start the simulation. The latter should be
-contained in all folders. 
+contained in all folders.
+
+### SOLVER CONFIG
+~~~python
+# GENERAL SOLVER CONFIG  ######################################################
+# maximal iteration per timestep
+max_iter_num = 500
+FEM_Lagrange_degree = 1
+~~~
+This block sets the maximal number of L-iterations `max_iter_num` that the
+LDD solver gets for the calculation of each time step.
+`FEM_Lagrange_degree` sets the degree of FEM ansatz functions that are used.
+While the communication of dofs over interfaces has been writen in a way to
+allow for general anasatz functions, other parts of the code have not.
+If higher order FEM methods are to be used, the code will have to be adjusted.
+
+### GRID AND MESH STUDY CONFIG
+~~~python
+# GRID AND MESH STUDY SPECIFICATIONS  #########################################
+mesh_study = True
+resolutions = {
+                # 1: 1e-6,
+                # 2: 1e-6,
+                # 4: 1e-6,
+                # 8: 1e-6,
+                # 16: 5e-6,
+                # 32: 5e-6,
+                64: 2e-6,
+                # 128: 1e-6,
+                # 256: 1e-6,
+                }
+
+# 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.0}
+timestep_size = 0.001
+number_of_timesteps = 20
+~~~
+This block sets mesh parameters, number of timesteps to be simulated, timestep
+size as well as whether or not a mesh study is being performed.
+In detail:
+- ~python mesh_study = True/False~ Toggle mesh study mode.
+  This changes the output slightly. Have a look at examples in the `mesh_study`
+  folders.
+- resolutions = {
+                # 1: 1e-6,
+                # 2: 1e-6,
+                # 4: 1e-6,
+                # 8: 1e-6,
+                # 16: 5e-6,
+                # 32: 5e-6,
+                64: 2e-6,
+                # 128: 1e-6,
+                # 256: 1e-6,
+                }
+
+- # 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.0}
+- timestep_size = 0.001
+- number_of_timesteps = 20
+
+### LDD SCHEME PARAMATERS
+~~~python
+# LDD scheme parameters  ######################################################
+Lw1 = 0.25 #/timestep_size
+Lnw1= 0.25
+
+Lw2 = 0.5 #/timestep_size
+Lnw2= 0.25
+
+lambda_w = 40
+lambda_nw = 40
+
+include_gravity = True
+debugflag = True
+analyse_condition = False
+~~~
+
+### I/O CONFIG
+~~~python
+# 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
+
+# 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,
+        # 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:
+    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
+    }
+~~~
+
+### OUTPUT FILE STRING
+~~~python
+# OUTPUT FILE STRING  #########################################################
+output_string = "./output/{}-{}_timesteps{}_P{}".format(
+    datestr, use_case, number_of_timesteps, FEM_Lagrange_degree
+    )
+~~~
+
+### DOMAIN SUBSTRUCTURING
+~~~python
+# DOMAIN AND INTERFACE  #######################################################
+substructuring = dss.twoSoilLayers()
+interface_def_points = substructuring.interface_def_points
+adjacent_subdomains = substructuring.adjacent_subdomains
+subdomain_def_points = substructuring.subdomain_def_points
+outer_boundary_def_points = substructuring.outer_boundary_def_points
+~~
+
+### MODEL CONFIGURATION
+~~~python
+# MODEL CONFIGURATION #########################################################
+isRichards = {
+    1: True, #
+    2: False
+    }
+
+
+viscosity = {#
+# subdom_num : viscosity
+    1: {'wetting' :1,
+         'nonwetting': 1/50}, #
+    2: {'wetting' :1,
+         'nonwetting': 1/50}
+}
+
+porosity = {#
+# subdom_num : porosity
+    1: 0.22,#
+    2: 0.22
+}
+
+# Dict of the form: { subdom_num : density }
+densities = {
+    1: {'wetting': 997,
+        'nonwetting': 1.225},
+    2: {'wetting': 997,
+        'nonwetting': 1.225}
+}
+
+gravity_acceleration = 9.81
+
+L = {#
+# subdom_num : subdomain L for L-scheme
+    1 : {'wetting' :Lw1,
+         'nonwetting': Lnw1},#
+    2 : {'wetting' :Lw2,
+         'nonwetting': Lnw2}
+}
+
+
+lambda_param = {#
+# subdom_num : lambda parameter for the L-scheme
+    0 : {'wetting' :lambda_w,
+         'nonwetting': lambda_nw},#
+}
+
+intrinsic_permeability = {
+    1: 0.01,
+    2: 0.01,
+}
+
+# relative permeabilties
+rel_perm_definition = {
+    1: {"wetting": "Spow2",
+        "nonwetting": "oneMinusSpow2"},
+    2: {"wetting": "Spow3",
+        "nonwetting": "oneMinusSpow3"},
+}
+
+rel_perm_dict = fts.generate_relative_permeability_dicts(rel_perm_definition)
+relative_permeability = rel_perm_dict["ka"]
+ka_prime = rel_perm_dict["ka_prime"]
+
+# S-pc relation
+Spc_on_subdomains = {
+    1: {"testSpc": {"index": 1}},
+    2: {"testSpc": {"index": 2}},
+}
+
+Spc = fts.generate_Spc_dicts(Spc_on_subdomains)
+S_pc_sym = Spc["symbolic"]
+S_pc_sym_prime = Spc["prime_symbolic"]
+sat_pressure_relationship = Spc["dolfin"]
+~~~
+
+### MANUFACTURED SOLUTION/SOURCE EXPRESSION
+~~~python
+###############################################################################
+# 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': (-7.0 - (1.0 + t*t)*(1.0 + x*x + y*y))}, #*(1-x)**2*(1+x)**2*(1-y)**2},
+    2: {'wetting': (-7.0 - (1.0 + t*t)*(1.0 + x*x)), #*(1-x)**2*(1+x)**2*(1+y)**2,
+        'nonwetting': (-2-t*(1.1+y + x**2))*y**2}, #*(1-x)**2*(1+x)**2*(1+y)**2},
+} #-y*y*(sym.sin(-2*t+2*x)*sym.sin(1/2*y-1.2*t)) - t*t*x*(0.5-y)*y*(1-x)
+~~~
+
+### BOUNDARY CONDITIONS
+~~~python
+# BOUNDARY CONDITIONS #########################################################
+# Dictionary of dirichlet boundary conditions. If an exact solution case is
+# used, use the hlp.generate_exact_DirichletBC() method to generate the
+# Dirichlet Boundary conditions from the exact solution.
+dirichletBC = hlp.generate_exact_DirichletBC(
+        isRichards=isRichards,
+        outer_boundary_def_points=outer_boundary_def_points,
+        exact_solution=exact_solution
+    )
+# If no exact solution is provided you need to provide a dictionary of boundary
+# conditions. See the definiton of hlp.generate_exact_DirichletBC() to see
+# the structure.
+~~~
diff --git a/Usecases/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-test.py b/Usecases/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-test.py
index 6836a27..5b2f496 100755
--- a/Usecases/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-test.py
+++ b/Usecases/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-test.py
@@ -201,177 +201,28 @@ intrinsic_permeability = {
     2: 0.01,
 }
 
-## 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
-
-_rel_perm1w = ft.partial(rel_perm1w)
-_rel_perm1nw = ft.partial(rel_perm1nw)
-
-subdomain1_rel_perm = {
-    'wetting': _rel_perm1w,#
-    'nonwetting': _rel_perm1nw
-}
-## relative permeabilty functions on subdomain 2
-def rel_perm2w(s):
-    # relative permeabilty wetting on subdomain2
-    return intrinsic_permeability[2]*s**3
-def rel_perm2nw(s):
-    # relative permeabilty nonwetting on subdomain2
-    return intrinsic_permeability[2]*(1-s)**3
-
-_rel_perm2w = ft.partial(rel_perm2w)
-_rel_perm2nw = ft.partial(rel_perm2nw)
-
-subdomain2_rel_perm = {
-    'wetting': _rel_perm2w,#
-    'nonwetting': _rel_perm2nw
-}
-
-## dictionary of relative permeabilties on all domains.
-relative_permeability = {#
-    1: subdomain1_rel_perm,
-    2: subdomain2_rel_perm
+# relative permeabilties
+rel_perm_definition = {
+    1: {"wetting": "Spow2",
+        "nonwetting": "oneMinusSpow2"},
+    2: {"wetting": "Spow3",
+        "nonwetting": "oneMinusSpow3"},
 }
 
+rel_perm_dict = fts.generate_relative_permeability_dicts(rel_perm_definition)
+relative_permeability = rel_perm_dict["ka"]
+ka_prime = rel_perm_dict["ka_prime"]
 
-# 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)
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 1
-def rel_perm2w_prime(s):
-    # relative permeabilty on subdomain2
-    return intrinsic_permeability[2]*3*s**2
-
-def rel_perm2nw_prime(s):
-    # relative permeabilty on subdomain2
-    return -3*intrinsic_permeability[2]*(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: subdomain2_rel_perm_prime,
-}
-
-
-# def saturation1(pc, subdomain_index):
-#     # inverse capillary pressure-saturation-relationship
-#     return df.conditional(pc > 0, 1/((1 + pc)**(1/(subdomain_index + 1))), 1)
-#
-# def saturation2(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 saturation1_sym(pc, subdomain_index):
-#     # inverse capillary pressure-saturation-relationship
-#     return 1/((1 + pc)**(1/(subdomain_index + 1)))
-#
-#
-# def saturation2_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 saturation1_sym_prime(pc, subdomain_index):
-#     # inverse capillary pressure-saturation-relationship
-#     return -(1/(subdomain_index + 1))*(1 + pc)**((-subdomain_index - 2)/(subdomain_index + 1))
-#
-#
-# def saturation2_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) )
-#
-# # 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(saturation1_sym, subdomain_index = 1),
-#     2: ft.partial(saturation2_sym, n_index=3, alpha=0.001),
-# }
-#
-# S_pc_sym_prime = {
-#     1: ft.partial(saturation1_sym_prime, subdomain_index = 1),
-#     2: ft.partial(saturation2_sym_prime, n_index=3, alpha=0.001),
-# }
-#
-# sat_pressure_relationship = {
-#     1: ft.partial(saturation1, subdomain_index = 1),#,
-#     2: ft.partial(saturation2, n_index=3, alpha=0.001),
-# }
-
-def saturation(pc, index):
-    # inverse capillary pressure-saturation-relationship
-    return df.conditional(pc > 0, 1/((1 + pc)**(1/(index + 1))), 1)
-
-
-def saturation_sym(pc, index):
-    # inverse capillary pressure-saturation-relationship
-    return 1/((1 + pc)**(1/(index + 1)))
-
-
-# derivative of S-pc relationship with respect to pc. This is needed for the
-# construction of a analytic solution.
-def saturation_sym_prime(pc, index):
-    # inverse capillary pressure-saturation-relationship
-    return -1/((index+1)*(1 + pc)**((index+2)/(index+1)))
-
-
-# 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, index=1),
-    2: ft.partial(saturation_sym, index=2),
-    # 3: ft.partial(saturation_sym, index=2),
-    # 4: ft.partial(saturation_sym, index=1)
-}
-
-S_pc_sym_prime = {
-    1: ft.partial(saturation_sym_prime, index=1),
-    2: ft.partial(saturation_sym_prime, index=2),
-    # 3: ft.partial(saturation_sym_prime, index=2),
-    # 4: ft.partial(saturation_sym_prime, index=1)
-}
-
-sat_pressure_relationship = {
-    1: ft.partial(saturation, index=1),
-    2: ft.partial(saturation, index=2),
-    # 3: ft.partial(saturation, index=2),
-    # 4: ft.partial(saturation, index=1)
+# S-pc relation
+Spc_on_subdomains = {
+    1: {"testSpc": {"index": 1}},
+    2: {"testSpc": {"index": 2}},
 }
 
+Spc = fts.generate_Spc_dicts(Spc_on_subdomains)
+S_pc_sym = Spc["symbolic"]
+S_pc_sym_prime = Spc["prime_symbolic"]
+sat_pressure_relationship = Spc["dolfin"]
 
 ###############################################################################
 # Manufacture source expressions with sympy #
-- 
GitLab