fix hardwired 0 in initial condition

LDDsimulation/helpers.py

@@ -72,12 +72,15 @@ def generate_exact_solution_expressions(
         dtpc = sym.diff(pc, t, 1)
         dxpc = sym.diff(pc, x, 1)
         dypc = sym.diff(pc, y, 1)
+        print(f" type of dtpc is {type(dtpc)}")
         if saturation_pressure_relationship is not None:
             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))
+            print(f" type of dS is {type(dS)}")
         else:
             S = symbolic_S_pc_relationship[subdomain]
             dS = symbolic_S_pc_relationship_prime[subdomain]
+            print(f" type of dS is {type(dS)}")
         for phase in subdomain_has_phases:
             # Turn above symbolic expression for exact solution into c code
             output['exact_solution'][subdomain].update(
@@ -85,7 +88,7 @@ def generate_exact_solution_expressions(
                 )
             # save the c code for initial conditions
             output['initial_condition'][subdomain].update(
-                {phase: sym.printing.ccode(symbolic_pressure[subdomain][phase].subs(t, 0))}
+                {phase: sym.printing.ccode(symbolic_pressure[subdomain][phase])}  #.subs(t, 0)
                 )
             if phase == "nonwetting":
                 dtS[subdomain].update(

TP-one-patch/debug_tests/R-one-patch-const-in-time.py

@@ -352,12 +352,6 @@ cutoff = gaussian/(gaussian + zero_on_shrinking)
 #         'nonwetting': (-1 -t*(1+y + x**2)**2)*cutoff},
 # }
 
-p_e_sym = {
-    0: {'wetting': -6 -(1 + x*x + y*y),
-        'nonwetting': -1 -(sym.sin(3*(1+y)/2*sym.pi)*sym.sin(5*(1+x)/2*sym.pi))**2},
-}
-
-
 print(f"\n\n\nsymbolic type is {type(p_e_sym[0]['wetting'])}\n\n\n")
 # # pw_sym_x*pw_sym_y
 # p_e_sym = {
@@ -370,6 +364,11 @@ print(f"\n\n\nsymbolic type is {type(p_e_sym[0]['wetting'])}\n\n\n")
 #         'nonwetting': -1*zero_on_shrinking+ 0*t},
 # }
 
+p_e_sym = {
+    0: {'wetting': -6 -(1 + x*x + y*y) + 0*t}
+        # 'nonwetting': -1 -(sym.sin(3*(1+y)/2*sym.pi)*sym.sin(5*(1+x)/2*sym.pi))**2},
+}
+
 
 pc_e_sym = dict()
 for subdomain, isR in isRichards.items():
@@ -396,11 +395,13 @@ for subdomain, isR in isRichards.items():
 #         )
 # }
 
+
 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(

TP-one-patch/mesh_study_for_fixed_timestep/TP-one-patch-mesh-study-fixed-timestep.py

@@ -34,13 +34,13 @@ resolutions = {
                 32: 5e-7,
                 64: 5e-7,
                 128: 5e-7,
-                # 256: 1e-10,
+                256: 5e-7,
                 # 512: 1e-10,
                 }
 
 ############ GRID #######################
 # mesh_resolution = 20
-timestep_size = 0.001
+timestep_size = 0.0025
 number_of_timesteps = 1
 plot_timestep_every = 1
 # decide how many timesteps you want analysed. Analysed means, that we write out