Merge branch 'master' of gitlab.mathematik.uni-stuttgart.de:seusdd/ldd-tpr-numerical-experiments

LDDsimulation/LDDsimulation.py

@@ -265,6 +265,7 @@ class LDDsimulation(object):
                 raise RuntimeError("include_gravity = True but no densities were given.")
         self._parameters_set = True
 
+
     def initialise(self) -> None:
         """ initialise LDD simulation
 

LDDsimulation/domainPatch.py

@@ -621,7 +621,7 @@ class DomainPatch(df.SubDomain):
                         # print(f"np.allclose({p_coord}, {gli_dof_coord}, rtol=1e-10, atol=1e-14): {b_close_a}")
                         if a_close_b and b_close_a:
                             p_previous_dof = p_prev_dof
-                    print("\n")
+                            
                     if p_previous_dof is None:
                         raise RuntimeError(f"didn't find p_previous_dof corresponding to dict key ({gli_dof_coord})",
                                             "something is wrong")

RR-2-patch-test-case/RR-2-patch-test.py

@@ -15,8 +15,8 @@ solver_tol = 1E-6
 
 ############ GRID #######################ü
 mesh_resolution = 40
-timestep_size = 0.01
-number_of_timesteps = 15
+timestep_size = 0.001
+number_of_timesteps = 1500
 # 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
@@ -25,11 +25,11 @@ starttime = 0
 Lw = 0.25 #/timestep_size
 Lnw=Lw
 
-l_param_w = 4
+lambda_w = 4
 
-include_gravity = False
+include_gravity = True
 debugflag = False
-analyse_condition = False
+analyse_condition = True
 
 output_string = "./output/new_gravity_term-number_of_timesteps{}_".format(number_of_timesteps)
 
@@ -138,8 +138,8 @@ L = {#
 
 lambda_param = {#
 # subdom_num : lambda parameter for the L-scheme
-    1 : {'wetting': l_param_w},#
-    2 : {'wetting': l_param_w}
+    1 : {'wetting': lambda_w},#
+    2 : {'wetting': lambda_w}
 }
 
 ## relative permeabilty functions on subdomain 1
@@ -226,7 +226,7 @@ write_to_file = {
 }
 
 # initialise LDD simulation class
-simulation = ldd.LDDsimulation(tol = 1E-14, debug = False, LDDsolver_tol=1E-9)
+simulation = ldd.LDDsimulation(tol = 1E-14, debug = debugflag, LDDsolver_tol=solver_tol)
 simulation.set_parameters(use_case=use_case,
     output_dir = output_string,#
     subdomain_def_points = subdomain_def_points,#
@@ -255,7 +255,7 @@ simulation.set_parameters(use_case=use_case,
     )
 
 simulation.initialise()
-print(simulation.__dict__)
+#print(simulation.__dict__)
 simulation.run(analyse_condition=analyse_condition)
 # simulation.LDDsolver(time = 0, debug = True, analyse_timestep = True)
 # df.info(parameters, True)

RR-multi-patch-with-inner-patch/RR-multi-patch-with-inner-patch.py

@@ -7,28 +7,34 @@ 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()
 
-# ----------------------------------------------------------------------------#
-# ------------------- MESH ---------------------------------------------------#
-# ----------------------------------------------------------------------------#
-mesh_resolution = 51
-# ----------------------------------------:-------------------------------------#
-# ------------------- TIME ---------------------------------------------------#
-# ----------------------------------------------------------------------------#
+use_case = "RR-multi-patch-with-inner-patch-cuttoff"
+solver_tol = 1E-6
+
+############ GRID #######################ü
+mesh_resolution = 20
 timestep_size = 0.01
 number_of_timesteps = 150
-# decide how many timesteps you want analysed. Analysed means, that we write
-# out subsequent errors of the L-iteration within the timestep.
+# 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 = 0.25/timestep_size
-Lw = 0.25
+Lw = 0.25 #/timestep_size
+Lnw=Lw
+
+lambda_w = 4
+
+include_gravity = True
+debugflag = False
+analyse_condition = False
 
+output_string = "./output/new_gravity_cutoff_term-number_of_timesteps{}_".format(number_of_timesteps)
 
 # ----------------------------------------------------------------------------#
 # ------------------- Domain and Interface -----------------------------------#
@@ -205,15 +211,15 @@ L = {
     5: {'wetting': Lw}
 }
 
-lamdal_w = 32
+lambda_w = 32
 
 # 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},
-    5: {'wetting': lamdal_w}
+    1: {'wetting': lambda_w},
+    2: {'wetting': lambda_w},
+    3: {'wetting': lambda_w},
+    4: {'wetting': lambda_w},
+    5: {'wetting': lambda_w}
 }
 
 
@@ -460,84 +466,37 @@ p_e_sym = {
 }
 
 
-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'],
-    5: -1*p_e_sym[5]['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']})
     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)}
+        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()
@@ -568,87 +527,15 @@ for subdomain in isRichards.keys():
                 )
 
 
-# # construction of the rhs that matches the above exact solution.
-# dtS = dict()
-# div_flux = dict()
-# source_expression = dict()
-# for subdomain, isR in isRichards.items():
-#     dtS.update({subdomain: dict()})
-#     div_flux.update({subdomain: dict()})
-#     source_expression.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:
-#         if phase == "nonwetting":
-#             dS = -dS
-#         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":
-#             dxdxflux = 1/mu*dka(1-S)*dS*dxpc*dxpa - 1/mu*dxdxpa*ka(1-S)
-#             dydyflux = 1/mu*dka(1-S)*dS*dypc*(dypa - rho*g) \
-#                 - 1/mu*dydypa*ka(1-S)
-#         else:
-#             dxdxflux = -1/mu*dka(S)*dS*dxpc*dxpa - 1/mu*dxdxpa*ka(S)
-#             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])
-
-# similarly to the outer boundary dictionary, if a patch has no outer boundary
-# None should be written instead of an expression.
-# This is a bit of a brainfuck:
-# dirichletBC[ind] gives a dictionary of the outer boundaries of subdomain ind.
-# Since a domain patch can have several disjoint outer boundary parts, the
-# expressions need to get an enumaration index which starts at 0.
-# So dirichletBC[ind][j] is the dictionary of outer dirichlet conditions of
-# subdomain ind and boundary part j.
-# Finally, dirichletBC[ind][j]['wetting'] and dirichletBC[ind][j]['nonwetting']
-# return the actual expression needed for the dirichlet condition for both
-# phases if present.
-
-# subdomain index: {outer boudary part index: {phase: expression}}
-# dirichletBC = {
-#     1: {0: {'wetting': exact_solution[1]['wetting']}},
-#     2: {0: {'wetting': exact_solution[2]['wetting']}},
-#     3: {0: {'wetting': exact_solution[3]['wetting']}},
-#     4: {0: {'wetting': exact_solution[4]['wetting']}},
-#     5: {0: {'wetting': exact_solution[5]['wetting']}}
-# }
-
 write_to_file = {
     'meshes_and_markers': True,
     'L_iterations': True
 }
 
 # initialise LDD simulation class
-simulation = ldd.LDDsimulation(tol=1E-14, debug=False, LDDsolver_tol=1E-7)
-simulation.set_parameters(output_dir="./output/with_cutoff_function",
+simulation = ldd.LDDsimulation(tol=1E-14, debug=debugflag, LDDsolver_tol=solver_tol)
+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,
@@ -670,12 +557,12 @@ simulation.set_parameters(output_dir="./output/with_cutoff_function",
                           dirichletBC_expression_strings=dirichletBC,
                           exact_solution=exact_solution,
                           densities=densities,
-                          include_gravity=True,
+                          include_gravity=include_gravity,
                           write2file=write_to_file,
                           )
 
 simulation.initialise()
 # print(simulation.__dict__)
-simulation.run()
+simulation.run(analyse_condition=analyse_condition)
 # simulation.LDDsolver(time=0, debug=True, analyse_timestep=True)
 # df.info(parameters, True)

TP-R-multi-patch-same-wetting-phase-as-RR-zero-nonwetting/TP-R-multi-patch-same-wetting-phase-as-RR-zero-nonwetting.py

@@ -17,7 +17,7 @@ use_case = "TP-R-multi-patch-zero-nonwetting-phase"
 solver_tol = 1E-6
 
 ############ GRID #######################ü
-mesh_resolution = 50
+mesh_resolution = 40
 timestep_size = 0.001
 number_of_timesteps = 1500
 # decide how many timesteps you want analysed. Analysed means, that we write out
@@ -25,7 +25,7 @@ number_of_timesteps = 1500
 number_of_timesteps_to_analyse = 10
 starttime = 0
 
-Lw = 0.25 #/timestep_size
+Lw = 0.5 #/timestep_size
 Lnw=Lw
 
 l_param_w = 40

TP-R-two-patch-test-case/TP-R-2-patch-test.py

@@ -13,29 +13,29 @@ import helpers as hlp
 # init sympy session
 sym.init_printing()
 
-use_case = "TP-R-two-patch"
-solver_tol = 5E-7
+use_case = "TP-R-two-patch-realistic    "
+solver_tol = 1E-6
 
 ############ GRID #######################ü
-mesh_resolution = 40
-timestep_size = 0.000001
+mesh_resolution = 20
+timestep_size = 0.00001
 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
+number_of_timesteps_to_analyse = 5
 starttime = 0
 
-Lw = 5 #/timestep_size
+Lw = 1 #/timestep_size
 Lnw=Lw
 
-l_param_w = 100
-l_param_nw = 100
+l_param_w = 40
+l_param_nw = 40
 
 include_gravity = True
 debugflag = True
 analyse_condition = False
 
-output_string = "./output/after_reimplementing_gravity_term_number_of_timesteps{}_".format(number_of_timesteps)
+output_string = "./output/after_reimplementing_gravity_term-realistic_number_of_timesteps{}_".format(number_of_timesteps)
 
 ##### Domain and Interface ####
 # global simulation domain domain

TP-one-patch/TP-one-patch-purely-postive-pc.py

@@ -17,10 +17,10 @@ sym.init_printing()
 solver_tol = 1E-6
 
 ############ GRID #######################
-
-mesh_resolution = 40
-timestep_size = 0.001
-number_of_timesteps = 10
+use_case="TP-one-patch-purely-positive-pc"
+mesh_resolution = 30
+timestep_size = 0.0005
+number_of_timesteps = 2500
 # 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 = 1