set up new TPR 2 dom mesh study for the new code. the previous ones were...

set up new TPR 2 dom mesh study for the new code. the previous ones were without the extra nonwetting terms

Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study.py

@@ -24,26 +24,26 @@ use_case = "TPR-2-patch-realistic"
 max_iter_num = 500
 FEM_Lagrange_degree = 1
 mesh_study = True
-resolutions = { 1: 7e-7,
-                2: 7e-7,
-                4: 7e-7,
-                8: 7e-7,
-                16: 7e-7,
-                32: 7e-7,
-                64: 7e-7,
-                128: 7e-7,
-                #256: 7e-7,
+resolutions = { 1: 1e-6,
+                2: 1e-6,
+                4: 1e-6,
+                8: 1e-6,
+                16: 1e-6,
+                32: 1e-6,
+                64: 1e-6,
+                128: 1e-6,
+                #256: 1e-6,
                 }
 
 ############ GRID #######################
 # mesh_resolution = 20
 timestep_size = 0.001
-number_of_timesteps = 800
+number_of_timesteps = 1000
 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 = 4
-starttimes = [0.0, 0.7]
+starttimes = [0.0]
 
 Lw = 0.025 #/timestep_size
 Lnw=Lw
@@ -152,7 +152,7 @@ isRichards = {
 viscosity = {#
 # subdom_num : viscosity
     1 : {'wetting' :1},
-         #'nonwetting': 1}, #
+         'nonwetting': 1/50}, #
     2 : {'wetting' :1,
          'nonwetting': 1/50}
 }
@@ -166,8 +166,9 @@ porosity = {#
 # Dict of the form: { subdom_num : density }
 densities = {
     1: {'wetting': 997},
-    2: {'wetting': 997,
         'nonwetting': 1.225},
+    2: {'wetting': 997,
+        'nonwetting': 1.225}
 }
 
 gravity_acceleration = 9.81
@@ -175,7 +176,7 @@ gravity_acceleration = 9.81
 L = {#
 # subdom_num : subdomain L for L-scheme
     1 : {'wetting' :Lw},
-         # 'nonwetting': 0.25},#
+         'nonwetting': Lnw},#
     2 : {'wetting' :Lw,
          'nonwetting': Lnw}
 }
@@ -189,28 +190,41 @@ lambda_param = {#
          'nonwetting': lambda_nw}
 }
 
+# intrinsic_permeability = {
+#     1: {"wetting": 1,
+#         "nonwetting": 1},
+#     2: {"wetting": 1,
+#         "nonwetting": 1},
+# }
+
+intrinsic_permeability = {
+    1: 1,
+    2: 1,
+}
+
 ## relative permeabilty functions on subdomain 1
 def rel_perm1w(s):
     # relative permeabilty wetting on subdomain1
-    return s**2
+    return intrinsic_permeability[1]*s**2
 
-# def rel_perm1nw(s):
-#     # relative permeabilty nonwetting on subdomain1
-#     return (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)
+_rel_perm1nw = ft.partial(rel_perm1nw)
+
 subdomain1_rel_perm = {
     'wetting': _rel_perm1w,#
-    # 'nonwetting': _rel_perm1nw
+    'nonwetting': _rel_perm1nw
 }
 ## relative permeabilty functions on subdomain 2
 def rel_perm2w(s):
     # relative permeabilty wetting on subdomain2
-    return s**3
+    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
-    return (1-s)**3
+    return intrinsic_permeability[2]*(1-s)**3
 
 _rel_perm2w = ft.partial(rel_perm2w)
 _rel_perm2nw = ft.partial(rel_perm2nw)
@@ -231,30 +245,30 @@ relative_permeability = {#
 # relative permeabilty functions on subdomain 1
 def rel_perm1w_prime(s):
     # relative permeabilty on subdomain1
-    return 2*s
+    return intrinsic_permeability[1]*2*s
 
-# def rel_perm1nw_prime(s):
-#     # relative permeabilty on subdomain1
-#     return 2*(1-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 subdomain1
-    return 3*s**2
+    # relative permeabilty on subdomain2
+    return intrinsic_permeability[2]*3*s**2
 
 def rel_perm2nw_prime(s):
-    # relative permeabilty on subdomain1
-    return -3*(1-s)**2
+    # 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_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
+    'wetting': _rel_perm1w_prime,
+    'nonwetting': _rel_perm1nw_prime
 }
 
 
@@ -441,6 +455,11 @@ for subdomain in isRichards.keys():
 #
 # sa
 
+f = open('TP-R-2-patch-mesh-study.py', 'r')
+print(f.read())
+f.close()
+
+
 for starttime in starttimes:
     for mesh_resolution, solver_tol in resolutions.items():
         # initialise LDD simulation class