set up TPR and TP-TP injection examples

Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-pure-dd.py

@@ -29,8 +29,8 @@ resolutions = {
                 # 2: 7e-7,
                 # 4: 7e-7,
                 # 8: 7e-7,
-                16: 5e-7,
-                # 32: 1e-7,
+                # 16: 5e-7,
+                32: 1e-7,
                 # 64: 7e-7,
                 # 128: 7e-7,
                 # 256: 7e-7
@@ -39,21 +39,22 @@ resolutions = {
 ############ GRID #######################
 # mesh_resolution = 20
 timestep_size = 0.005
-number_of_timesteps = 25
+number_of_timesteps = 200
+
 plot_timestep_every = 1
 # 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 = 0
 starttimes = [0.0]
 
-Lw = 0.025 #/timestep_size
-Lnw= 0.025
+Lw = 0.05 #/timestep_size
+Lnw= 0.05
 
 lambda_w = 20
 lambda_nw = 20
 include_gravity = True
 debugflag = False
-analyse_condition = False
+analyse_condition = True
 
 output_string = "./output/{}-{}_timesteps{}_P{}".format(datestr, use_case, number_of_timesteps, FEM_Lagrange_degree)
 

Two-phase-Two-phase/multi-patch/TP-TP-layered-soil-case-with-inner-patch/TP-TP-layered_soil_with_inner_patch-realistic-pure-dd.py

@@ -44,23 +44,23 @@ resolutions = {
 
 ############ GRID #######################
 # mesh_resolution = 20
-timestep_size = 0.00001
+timestep_size = 0.005
 number_of_timesteps = 20
-plot_timestep_every = 5
+plot_timestep_every = 1
 # 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 = 1
 starttime = 0.0
 
-Lw = 0.5 #/timestep_size
+Lw = 0.05 #/timestep_size
 Lnw=Lw
 
-lambda_w = 400
-lambda_nw = 400
+lambda_w = 40
+lambda_nw = 40
 
 include_gravity = True
 debugflag = True
-analyse_condition = True
+analyse_condition = False
 
 if mesh_study:
     output_string = "./output/{}-{}_timesteps{}_P{}".format(datestr, use_case, number_of_timesteps, FEM_Lagrange_degree)

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

@@ -13,31 +13,80 @@ import helpers as hlp
 # init sympy session
 sym.init_printing()
 
+date = datetime.datetime.now()
+datestr = date.strftime("%Y-%m-%d")
+#import ufl as ufl
 
-solver_tol = 4E-6
+# init sympy session
+sym.init_printing()
 
-############ GRID #######################
 use_case = "TP-one-patch-purely-positive-pc"
-mesh_resolution = 60
+# solver_tol = 5E-7
+max_iter_num = 1000
+FEM_Lagrange_degree = 1
+mesh_study = False
+resolutions = {
+                # 1: 1e-7,  # h=2
+                # 2: 2e-5,  # h=1.1180
+                # 4: 1e-6,  # h=0.5590
+                # 8: 1e-6,  # h=0.2814
+                # 16: 5e-7, # h=0.1412
+                # 32: 1e-6,
+                64: 5e-7,
+                # 128: 5e-7
+                }
+
+
+############ GRID #######################
+# mesh_resolution = 20
 timestep_size = 0.0005
 number_of_timesteps = 2500
+plot_timestep_every = 1
 # 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 = 5
-starttime = 0
+number_of_timesteps_to_analyse = 0
+starttime = 0.0
 
 Lw = 0.25 #/timestep_size
 Lnw=Lw
 
-l_param_w = 40
-l_param_nw = 40
+lambda_w = 40
+lambda_nw = 40
 
-include_gravity = True
-debugflag = False
-analyse_condition = True
+include_gravity = False
+debugflag = True
+analyse_condition = False
+
+if mesh_study:
+    output_string = "./output/{}-{}_timesteps{}_P{}".format(datestr, use_case, number_of_timesteps, FEM_Lagrange_degree)
+else:
+    for tol in resolutions.values():
+        solver_tol = tol
+    output_string = "./output/{}-{}_timesteps{}_P{}_solver_tol{}".format(datestr, use_case, number_of_timesteps, FEM_Lagrange_degree, solver_tol)
+
+
+# toggle what should be written to files
+if mesh_study:
+    write_to_file = {
+        'space_errornorms': True,
+        'meshes_and_markers': True,
+        'L_iterations_per_timestep': False,
+        'solutions': False,
+        'absolute_differences': False,
+        'condition_numbers': analyse_condition,
+        'subsequent_errors': False
+    }
+else:
+    write_to_file = {
+        'space_errornorms': True,
+        'meshes_and_markers': True,
+        'L_iterations_per_timestep': True,
+        'solutions': True,
+        'absolute_differences': True,
+        'condition_numbers': analyse_condition,
+        'subsequent_errors': True
+    }
 
-use_case="TP-one-patch-purely-positive-pc"
-output_string = "./output/2019-08-26-pure_positive_pc_number_of_timesteps{}_".format(number_of_timesteps)
 
 ##### Domain and Interface ####
 # global simulation domain domain
@@ -106,8 +155,8 @@ L = {#
 
 lambda_param = {#
 # subdom_num : lambda parameter for the L-scheme
-    0: {'wetting' :l_param_w,
-         'nonwetting': l_param_nw},#
+    0: {'wetting' :lambda_w,
+         'nonwetting': lambda_nw},#
 }
 
 ## relative permeabilty functions on subdomain 1
@@ -373,46 +422,72 @@ for subdomain in isRichards.keys():
                 )
 
 
-# def saturation(pressure, subdomain_index):
-#     # inverse capillary pressure-saturation-relationship
-#     return df.conditional(pressure < 0, 1/((1 - pressure)**(1/(subdomain_index + 1))), 1)
-#
-# sa
-
-write_to_file = {
-    'meshes_and_markers': True,
-    'L_iterations': True
-}
 
+for mesh_resolution, solver_tol in resolutions.items():
     # initialise LDD simulation class
-simulation = ldd.LDDsimulation(tol = 1E-14, LDDsolver_tol=solver_tol, debug=debugflag)
+    simulation = ldd.LDDsimulation(
+        tol=1E-14,
+        LDDsolver_tol=solver_tol,
+        debug=debugflag,
+        max_iter_num=max_iter_num,
+        FEM_Lagrange_degree=FEM_Lagrange_degree,
+        mesh_study=mesh_study
+        )
+
     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,#
-    outer_boundary_def_points = outer_boundary_def_points,#
-    adjacent_subdomains = adjacent_subdomains,#
-    mesh_resolution = mesh_resolution,#
-    viscosity = viscosity,#
-    porosity = porosity,#
-    L = L,#
-    lambda_param = lambda_param,#
-    relative_permeability = relative_permeability,#
-    saturation = sat_pressure_relationship,#
-    starttime = starttime,#
+                              output_dir=output_string,
+                              subdomain_def_points=subdomain_def_points,
+                              isRichards=isRichards,
+                              interface_def_points=interface_def_points,
+                              outer_boundary_def_points=outer_boundary_def_points,
+                              adjacent_subdomains=adjacent_subdomains,
+                              mesh_resolution=mesh_resolution,
+                              viscosity=viscosity,
+                              porosity=porosity,
+                              L=L,
+                              lambda_param=lambda_param,
+                              relative_permeability=relative_permeability,
+                              saturation=sat_pressure_relationship,
+                              starttime=starttime,
                               number_of_timesteps=number_of_timesteps,
                               number_of_timesteps_to_analyse=number_of_timesteps_to_analyse,
-    timestep_size = timestep_size,#
-    sources = source_expression,#
-    initial_conditions = initial_condition,#
-    dirichletBC_expression_strings = dirichletBC,#
-    exact_solution = exact_solution,#
+                              plot_timestep_every=plot_timestep_every,
+                              timestep_size=timestep_size,
+                              sources=source_expression,
+                              initial_conditions=initial_condition,
+                              dirichletBC_expression_strings=dirichletBC,
+                              exact_solution=exact_solution,
                               densities=densities,
                               include_gravity=include_gravity,
-    write2file = write_to_file,#
+                              write2file=write_to_file,
                               )
 
     simulation.initialise()
+    output_dir = simulation.output_dir
     # simulation.write_exact_solution_to_xdmf()
-simulation.run(analyse_condition=analyse_condition)
+    output = simulation.run(analyse_condition=analyse_condition)
+    for subdomain_index, subdomain_output in output.items():
+        mesh_h = subdomain_output['mesh_size']
+        for phase, different_errornorms in subdomain_output['errornorm'].items():
+            filename = output_dir + "subdomain{}-space-time-errornorm-{}-phase.csv".format(subdomain_index, phase)
+            # for errortype, errornorm in different_errornorms.items():
+
+                # eocfile = open("eoc_filename", "a")
+                # eocfile.write( str(mesh_h) + " " + str(errornorm) + "\n" )
+                # eocfile.close()
+                # if subdomain.isRichards:mesh_h
+            data_dict = {
+                'mesh_parameter': mesh_resolution,
+                'mesh_h': mesh_h,
+            }
+            for error_type, errornorms in different_errornorms.items():
+                data_dict.update(
+                    {error_type: errornorms}
+                )
+            errors = pd.DataFrame(data_dict, index=[mesh_resolution])
+            # check if file exists
+            if os.path.isfile(filename) == True:
+                with open(filename, 'a') as f:
+                    errors.to_csv(f, header=False, sep='\t', encoding='utf-8', index=False)
+            else:
+                errors.to_csv(filename, sep='\t', encoding='utf-8', index=False)