fix weird git fuckug

RR-multi-patch-plus-gravity-const-solution/RR-multi-patch-with-gravity-constant-solution.py

@@ -7,25 +7,33 @@ import sympy as sym
 # 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()
 
+solver_tol = 1E-7
+
 # ----------------------------------------------------------------------------#
 # ------------------- MESH ---------------------------------------------------#
 # ----------------------------------------------------------------------------#
-mesh_resolution = 40
+mesh_resolution = 10
 # ----------------------------------------:-------------------------------------#
 # ------------------- TIME ---------------------------------------------------#
 # ----------------------------------------------------------------------------#
 timestep_size = 0.01
-number_of_timesteps = 100
+number_of_timesteps = 50
 # 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
+number_of_timesteps_to_analyse = 10
 starttime = 0
 
+include_gravity = True
+Lw = 0.25
+lambda_w = 4
+
+output_string = "./output/testing_number_of_timesteps_{}_".format(number_of_timesteps)
 
 # ----------------------------------------------------------------------------#
 # ------------------- Domain and Interface -----------------------------------#
@@ -153,6 +161,7 @@ densities = {
 }
 
 gravity_acceleration = 9.81
+
 # Dict of the form: { subdom_num : porosity }
 porosity = {
     1: 1,
@@ -163,19 +172,19 @@ porosity = {
 
 # subdom_num : subdomain L for L-scheme
 L = {
-    1: {'wetting': 0.25},
-    2: {'wetting': 0.25},
-    3: {'wetting': 0.25},
-    4: {'wetting': 0.25}
+    1: {'wetting': Lw},
+    2: {'wetting': Lw},
+    3: {'wetting': Lw},
+    4: {'wetting': Lw}
 }
 
-lamdal_w = 4
+
 # subdom_num : lambda parameter for the L-scheme
 lambda_param = {
-    1: {'wetting': lamdal_w},
-    2: {'wetting': lamdal_w},
-    3: {'wetting': lamdal_w},
-    4: {'wetting': lamdal_w}
+    1: {'wetting': lambda_w},
+    2: {'wetting': lambda_w},
+    3: {'wetting': lambda_w},
+    4: {'wetting': lambda_w}
 }
 
 
@@ -321,82 +330,47 @@ p_e_sym = {
     4: {'wetting': -3 + 0*t}
 }
 
-pc_e_sym = {
-    1: -1*p_e_sym[1]['wetting'],
-    2: -1*p_e_sym[2]['wetting'],
-    3: -1*p_e_sym[3]['wetting'],
-    4: -1*p_e_sym[4]['wetting']
-}
+# pc_e_sym = {
+#     1: -1*p_e_sym[1]['wetting'],
+#     2: -1*p_e_sym[2]['wetting'],
+#     3: -1*p_e_sym[3]['wetting'],
+#     4: -1*p_e_sym[4]['wetting']
+# }
 
-# construction of 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"]
-    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)
+        pc_e_sym.update({subdomain: -p_e_sym[subdomain]['wetting']})
     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)}
+        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,
                         )
-        # print(f"source_expression[{subdomain}][{phase}] =", source_expression[subdomain][phase])
+
+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()
@@ -434,8 +408,8 @@ write_to_file = {
 }
 
 # initialise LDD simulation class
-simulation = ldd.LDDsimulation(tol=1E-14, debug=True, LDDsolver_tol=5E-7)
-simulation.set_parameters(output_dir="./output/",
+simulation = ldd.LDDsimulation(tol=1E-14, debug=False, LDDsolver_tol=solver_tol)
+simulation.set_parameters(output_dir=output_string,
                           subdomain_def_points=subdomain_def_points,
                           isRichards=isRichards,
                           interface_def_points=interface_def_points,
@@ -457,7 +431,7 @@ 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,
                           )