From 60d6185f65e09bac3358004ac0f7fb97c72fc5a7 Mon Sep 17 00:00:00 2001
From: David Seus <david.seus@ians.uni-stuttgart.de>
Date: Fri, 16 Aug 2019 15:07:14 +0200
Subject: [PATCH] fix weird git fuckug

---
 TP-R-two-patch-test-case/TP-R-2-patch-test.py | 178 ++++++++----------
 1 file changed, 76 insertions(+), 102 deletions(-)

diff --git a/TP-R-two-patch-test-case/TP-R-2-patch-test.py b/TP-R-two-patch-test-case/TP-R-2-patch-test.py
index 459f2dd..168870d 100755
--- a/TP-R-two-patch-test-case/TP-R-2-patch-test.py
+++ b/TP-R-two-patch-test-case/TP-R-2-patch-test.py
@@ -7,11 +7,36 @@ import typing as tp
 import domainPatch as dp
 import LDDsimulation as ldd
 import functools as ft
+import helpers as hlp
 #import ufl as ufl
 
 # init sympy session
 sym.init_printing()
 
+use_case = "TP-R-two-patch"
+solver_tol = 5E-7
+
+############ GRID #######################ü
+mesh_resolution = 40
+timestep_size = 0.000001
+number_of_timesteps = 15
+# 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
+
+Lw = 5 #/timestep_size
+Lnw=Lw
+
+l_param_w = 100
+l_param_nw = 100
+
+include_gravity = True
+debugflag = True
+analyse_condition = False
+
+output_string = "./output/after_reimplementing_gravity_term_number_of_timesteps{}_".format(number_of_timesteps)
+
 ##### Domain and Interface ####
 # global simulation domain domain
 sub_domain0_vertices = [df.Point(-1.0, -1.0),
@@ -80,53 +105,44 @@ isRichards = {
     }
 
 
-############ GRID #######################ü
-mesh_resolution = 50
-timestep_size = 0.01
-number_of_timesteps = 160
-# 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 = 11
-starttime = 0
-
 viscosity = {#
 # subdom_num : viscosity
     1 : {'wetting' :1},
          #'nonwetting': 1}, #
     2 : {'wetting' :1,
-         'nonwetting': 1}
+         'nonwetting': 1/50}
 }
 
 porosity = {#
 # subdom_num : porosity
-    1 : 1,#0.22,#
-    2 : 1#0.022
+    1 : 0.22,#
+    2 : 0.0022
 }
 
 # Dict of the form: { subdom_num : density }
 densities = {
-    1: {'wetting': 1},
-    2: {'wetting': 1,
-        'nonwetting': 1},
+    1: {'wetting': 997},
+    2: {'wetting': 997,
+        'nonwetting': 1.225},
 }
 
 gravity_acceleration = 9.81
 
 L = {#
 # subdom_num : subdomain L for L-scheme
-    1 : {'wetting' :0.25},
+    1 : {'wetting' :Lw},
          # 'nonwetting': 0.25},#
-    2 : {'wetting' :0.25,
-         'nonwetting': 0.25}
+    2 : {'wetting' :Lw,
+         'nonwetting': Lnw}
 }
 
-l_param = 40
+
 lambda_param = {#
 # subdom_num : lambda parameter for the L-scheme
-    1 : {'wetting' :l_param},
+    1 : {'wetting' :l_param_w},
          # 'nonwetting': l_param},#
-    2 : {'wetting' :l_param,
-         'nonwetting': l_param}
+    2 : {'wetting' :l_param_w,
+         'nonwetting': l_param_nw}
 }
 
 ## relative permeabilty functions on subdomain 1
@@ -185,7 +201,7 @@ def rel_perm2w_prime(s):
 
 def rel_perm2nw_prime(s):
     # relative permeabilty on subdomain1
-    return 3*(1-s)**2
+    return -3*(1-s)**2
 
 _rel_perm1w_prime = ft.partial(rel_perm1w_prime)
 # _rel_perm1nw_prime = ft.partial(rel_perm1nw_prime)
@@ -307,87 +323,44 @@ x, y = sym.symbols('x[0], x[1]')  # needed by UFL
 t = sym.symbols('t', positive=True)
 
 p_e_sym = {
-    1: {'wetting': 1.0 - (1.0 + t*t)*(1.0 + x*x + y*y)},
-    2: {'wetting': 1.0 - (1.0 + t*t)*(1.0 + x*x),
-        'nonwetting': (-t*(1-y - x**2)**2 - sym.sqrt(2+t**2))*y},
+    1: {'wetting': (-5.0 - (1.0 + t*t)*(1.0 + x*x + y*y))}, #*(1-x)**2*(1+x)**2*(1-y)**2},
+    2: {'wetting': (-5.0 - (1.0 + t*t)*(1.0 + x*x)), #*(1-x)**2*(1+x)**2*(1+y)**2,
+        'nonwetting': (-1-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)
 
 
-pc_e_sym = {
-    1: -1*p_e_sym[1]['wetting'],
-    2: p_e_sym[2]['nonwetting'] - p_e_sym[2]['wetting']
-}
-
-# turn above symbolic code into exact solution for dolphin and
-# construct the rhs that matches the above exact solution.
-dtS = dict()
-div_flux = dict()
-source_expression = dict()
-exact_solution = dict()
-initial_condition = dict()
+pc_e_sym = dict()
 for subdomain, isR in isRichards.items():
-    dtS.update({subdomain: dict()})
-    div_flux.update({subdomain: dict()})
-    source_expression.update({subdomain: dict()})
-    exact_solution.update({subdomain: dict()})
-    initial_condition.update({subdomain: dict()})
     if isR:
-        subdomain_has_phases = ["wetting"]
+        pc_e_sym.update({subdomain: -p_e_sym[subdomain]['wetting'].copy()})
     else:
-        subdomain_has_phases = ["wetting", "nonwetting"]
-
-    # conditional for S_pc_prime
-    pc = pc_e_sym[subdomain]
-    dtpc = sym.diff(pc, t, 1)
-    dxpc = sym.diff(pc, x, 1)
-    dypc = sym.diff(pc, y, 1)
-    S = sym.Piecewise((S_pc_sym[subdomain](pc), pc > 0), (1, True))
-    dS = sym.Piecewise((S_pc_sym_prime[subdomain](pc), pc > 0), (0, True))
-    for phase in subdomain_has_phases:
-        # Turn above symbolic expression for exact solution into c code
-        exact_solution[subdomain].update(
-            {phase: sym.printing.ccode(p_e_sym[subdomain][phase])}
-            )
-        # save the c code for initial conditions
-        initial_condition[subdomain].update(
-            {phase: sym.printing.ccode(p_e_sym[subdomain][phase].subs(t, 0))}
-            )
-        if phase == "nonwetting":
-            dtS[subdomain].update(
-                {phase: -porosity[subdomain]*dS*dtpc}
-                )
-        else:
-            dtS[subdomain].update(
-                {phase: porosity[subdomain]*dS*dtpc}
-                )
-        pa = p_e_sym[subdomain][phase]
-        dxpa = sym.diff(pa, x, 1)
-        dxdxpa = sym.diff(pa, x, 2)
-        dypa = sym.diff(pa, y, 1)
-        dydypa = sym.diff(pa, y, 2)
-        mu = viscosity[subdomain][phase]
-        ka = relative_permeability[subdomain][phase]
-        dka = ka_prime[subdomain][phase]
-        rho = densities[subdomain][phase]
-        g = gravity_acceleration
-
-        if phase == "nonwetting":
-            # x part of div(flux) for nonwetting
-            dxdxflux = -1/mu*dka(1-S)*dS*dxpc*dxpa + 1/mu*dxdxpa*ka(1-S)
-            # y part of div(flux) for nonwetting
-            dydyflux = -1/mu*dka(1-S)*dS*dypc*(dypa - rho*g) \
-                + 1/mu*dydypa*ka(1-S)
-        else:
-            # x part of div(flux) for wetting
-            dxdxflux = 1/mu*dka(S)*dS*dxpc*dxpa + 1/mu*dxdxpa*ka(S)
-            # y part of div(flux) for wetting
-            dydyflux = 1/mu*dka(S)*dS*dypc*(dypa - rho*g) + 1/mu*dydypa*ka(S)
-        div_flux[subdomain].update({phase: dxdxflux + dydyflux})
-        contructed_rhs = dtS[subdomain][phase] - div_flux[subdomain][phase]
-        source_expression[subdomain].update(
-            {phase: sym.printing.ccode(contructed_rhs)}
-            )
-        # print(f"source_expression[{subdomain}][{phase}] =", source_expression[subdomain][phase])
+        pc_e_sym.update({subdomain: p_e_sym[subdomain]['nonwetting'].copy()
+                                        - p_e_sym[subdomain]['wetting'].copy()})
+
+
+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()
@@ -431,8 +404,9 @@ write_to_file = {
 
 
 # initialise LDD simulation class
-simulation = ldd.LDDsimulation(tol = 1E-14, LDDsolver_tol = 1E-7, debug = False)
-simulation.set_parameters(output_dir = "./output/",#
+simulation = ldd.LDDsimulation(tol = 1E-14, LDDsolver_tol=solver_tol, debug=debugflag)
+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,#
@@ -454,10 +428,10 @@ simulation.set_parameters(output_dir = "./output/",#
     dirichletBC_expression_strings = dirichletBC,#
     exact_solution = exact_solution,#
     densities=densities,
-    include_gravity=True,
+    include_gravity=include_gravity,
     write2file = write_to_file,#
     )
 
 simulation.initialise()
 # simulation.write_exact_solution_to_xdmf()
-simulation.run()
+simulation.run(analyse_condition=analyse_condition)
-- 
GitLab