prepare tp-tp with garvity examples

Two-phase-Two-phase/two-patch/TP-TP-2-patch-pure-dd-avoid-interface-at-origin/TP-TP-2-patch-pure-dd-horizontal-interface-avoiding-origin.py

@@ -9,6 +9,8 @@ import LDDsimulation as ldd
 import functools as ft
 import helpers as hlp
 import datetime
+import os
+import pandas as pd
 
 date = datetime.datetime.now()
 datestr = date.strftime("%Y-%m-%d")
@@ -18,30 +20,71 @@ datestr = date.strftime("%Y-%m-%d")
 sym.init_printing()
 
 use_case = "TP-TP-2-patch-really-pure-dd-horizontal-interface-avoding-origin"
-solver_tol = 1E-6
+# 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 = 10
-timestep_size = 0.0001
-number_of_timesteps = 4000
+# mesh_resolution = 20
+timestep_size = 0.0005
+number_of_timesteps = 10
+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 = 6
-starttime = 0
+number_of_timesteps_to_analyse = 1
+starttime = 0.0
 
-Lw = 0.25 #/timestep_size
+Lw = 0.05 #/timestep_size
 Lnw=Lw
 
-lambda_w = 4
-lambda_nw = 4
+lambda_w = 400
+lambda_nw = 400
 
-include_gravity = False
+include_gravity = True
 debugflag = False
 analyse_condition = True
 
-output_string = "./output/{}-{}_timesteps{}_".format(datestr, use_case, number_of_timesteps)
+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': False,
+        'solutions': True,
+        'absolute_differences': True,
+        'condition_numbers': analyse_condition,
+        'subsequent_errors': True
+    }
 
 ##### Domain and Interface ####
 # global simulation domain domain
@@ -503,25 +546,15 @@ 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,
         max_iter_num=max_iter_num,
-    FEM_Lagrange_degree=FEM_Lagrange_degree
+        FEM_Lagrange_degree=FEM_Lagrange_degree,
+        mesh_study=mesh_study
         )
 
     simulation.set_parameters(use_case=use_case,
@@ -541,6 +574,7 @@ simulation.set_parameters(use_case=use_case,
                               starttime=starttime,
                               number_of_timesteps=number_of_timesteps,
                               number_of_timesteps_to_analyse=number_of_timesteps_to_analyse,
+                              plot_timestep_every=plot_timestep_every,
                               timestep_size=timestep_size,
                               sources=source_expression,
                               initial_conditions=initial_condition,
@@ -552,5 +586,31 @@ simulation.set_parameters(use_case=use_case,
                               )
 
     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)

Two-phase-Two-phase/two-patch/TP-TP-2-patch-test-case/TP-TP-2-patch-test.py

@@ -19,7 +19,7 @@ datestr = date.strftime("%Y-%m-%d")
 # init sympy session
 sym.init_printing()
 
-use_case = "TP-TP-2-patch-big-timestep"
+use_case = "TP-TP-2-patch"
 # solver_tol = 5E-7
 max_iter_num = 1000
 FEM_Lagrange_degree = 1
@@ -39,12 +39,12 @@ resolutions = {
 ############ GRID #######################
 # mesh_resolution = 20
 timestep_size = 0.001
-number_of_timesteps = 1
+number_of_timesteps = 10
 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 = 1
-starttime = 0.5
+starttime = 0.0
 
 Lw = 0.03 #/timestep_size
 Lnw=Lw