diff --git a/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 b/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
index e3a5a8e9c079712c14aa4ea75d958350e91fbfdd..25ce74e7e80180f6794eac61e56ac71d7aba9bc2 100755
--- a/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
+++ b/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,54 +546,71 @@ 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
-}
-
-
-# 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
- )
-
-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,
- 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,
- densities=densities,
- include_gravity=include_gravity,
- write2file=write_to_file,
- )
-
-simulation.initialise()
-# simulation.write_exact_solution_to_xdmf()
-simulation.run(analyse_condition=analyse_condition)
+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,
+ 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,
+ 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,
+ dirichletBC_expression_strings=dirichletBC,
+ exact_solution=exact_solution,
+ densities=densities,
+ include_gravity=include_gravity,
+ write2file=write_to_file,
+ )
+
+ simulation.initialise()
+ output_dir = simulation.output_dir
+ # simulation.write_exact_solution_to_xdmf()
+ 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)
diff --git a/Two-phase-Two-phase/two-patch/TP-TP-2-patch-test-case/TP-TP-2-patch-test.py b/Two-phase-Two-phase/two-patch/TP-TP-2-patch-test-case/TP-TP-2-patch-test.py
index 15937d60d0a5e0204bb3e11ff8089bc73cfb613d..31f42519805880f3354041a3b5d7efbbeeccd9eb 100755
--- a/Two-phase-Two-phase/two-patch/TP-TP-2-patch-test-case/TP-TP-2-patch-test.py
+++ b/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