diff --git a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study.py index 0517c79e84f2734501b759c6f6b8260efd04fb68..f17fddb00d9dca0a5efcf704850800eef45fcf56 100755 --- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study.py +++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study.py @@ -24,26 +24,26 @@ use_case = "TPR-2-patch-realistic" max_iter_num = 500 FEM_Lagrange_degree = 1 mesh_study = True -resolutions = { 1: 7e-7, - 2: 7e-7, - 4: 7e-7, - 8: 7e-7, - 16: 7e-7, - 32: 7e-7, - 64: 7e-7, - 128: 7e-7, - #256: 7e-7, +resolutions = { 1: 1e-6, + 2: 1e-6, + 4: 1e-6, + 8: 1e-6, + 16: 1e-6, + 32: 1e-6, + 64: 1e-6, + 128: 1e-6, + #256: 1e-6, } ############ GRID ####################### # mesh_resolution = 20 timestep_size = 0.001 -number_of_timesteps = 800 +number_of_timesteps = 1000 plot_timestep_every = 2 # 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 = 4 -starttimes = [0.0, 0.7] +starttimes = [0.0] Lw = 0.025 #/timestep_size Lnw=Lw @@ -152,7 +152,7 @@ isRichards = { viscosity = {# # subdom_num : viscosity 1 : {'wetting' :1}, - #'nonwetting': 1}, # + 'nonwetting': 1/50}, # 2 : {'wetting' :1, 'nonwetting': 1/50} } @@ -166,8 +166,9 @@ porosity = {# # Dict of the form: { subdom_num : density } densities = { 1: {'wetting': 997}, - 2: {'wetting': 997, 'nonwetting': 1.225}, + 2: {'wetting': 997, + 'nonwetting': 1.225} } gravity_acceleration = 9.81 @@ -175,7 +176,7 @@ gravity_acceleration = 9.81 L = {# # subdom_num : subdomain L for L-scheme 1 : {'wetting' :Lw}, - # 'nonwetting': 0.25},# + 'nonwetting': Lnw},# 2 : {'wetting' :Lw, 'nonwetting': Lnw} } @@ -189,28 +190,41 @@ lambda_param = {# 'nonwetting': lambda_nw} } +# intrinsic_permeability = { +# 1: {"wetting": 1, +# "nonwetting": 1}, +# 2: {"wetting": 1, +# "nonwetting": 1}, +# } + +intrinsic_permeability = { + 1: 1, + 2: 1, +} + ## relative permeabilty functions on subdomain 1 def rel_perm1w(s): # relative permeabilty wetting on subdomain1 - return s**2 + return intrinsic_permeability[1]*s**2 -# def rel_perm1nw(s): -# # relative permeabilty nonwetting on subdomain1 -# return (1-s)**2 +def rel_perm1nw(s): + # relative permeabilty nonwetting on subdomain1 + return intrinsic_permeability[1]*(1-s)**2 _rel_perm1w = ft.partial(rel_perm1w) -# _rel_perm1nw = ft.partial(rel_perm1nw) +_rel_perm1nw = ft.partial(rel_perm1nw) + subdomain1_rel_perm = { 'wetting': _rel_perm1w,# - # 'nonwetting': _rel_perm1nw + 'nonwetting': _rel_perm1nw } ## relative permeabilty functions on subdomain 2 def rel_perm2w(s): # relative permeabilty wetting on subdomain2 - return s**3 + return intrinsic_permeability[2]*s**3 def rel_perm2nw(s): # relative permeabilty nonwetting on subdosym.cos(0.8*t - (0.8*x + 1/7*y))main2 - return (1-s)**3 + return intrinsic_permeability[2]*(1-s)**3 _rel_perm2w = ft.partial(rel_perm2w) _rel_perm2nw = ft.partial(rel_perm2nw) @@ -231,30 +245,30 @@ relative_permeability = {# # relative permeabilty functions on subdomain 1 def rel_perm1w_prime(s): # relative permeabilty on subdomain1 - return 2*s + return intrinsic_permeability[1]*2*s -# def rel_perm1nw_prime(s): -# # relative permeabilty on subdomain1 -# return 2*(1-s) +def rel_perm1nw_prime(s): + # relative permeabilty on subdomain1 + return -1*intrinsic_permeability[1]*2*(1-s) # definition of the derivatives of the relative permeabilities # relative permeabilty functions on subdomain 1 def rel_perm2w_prime(s): - # relative permeabilty on subdomain1 - return 3*s**2 + # relative permeabilty on subdomain2 + return intrinsic_permeability[2]*3*s**2 def rel_perm2nw_prime(s): - # relative permeabilty on subdomain1 - return -3*(1-s)**2 + # relative permeabilty on subdomain2 + return -3*intrinsic_permeability[2]*(1-s)**2 _rel_perm1w_prime = ft.partial(rel_perm1w_prime) -# _rel_perm1nw_prime = ft.partial(rel_perm1nw_prime) +_rel_perm1nw_prime = ft.partial(rel_perm1nw_prime) _rel_perm2w_prime = ft.partial(rel_perm2w_prime) _rel_perm2nw_prime = ft.partial(rel_perm2nw_prime) subdomain1_rel_perm_prime = { - 'wetting': _rel_perm1w_prime - # 'nonwetting': _rel_perm1nw_prime + 'wetting': _rel_perm1w_prime, + 'nonwetting': _rel_perm1nw_prime } @@ -441,6 +455,11 @@ for subdomain in isRichards.keys(): # # sa +f = open('TP-R-2-patch-mesh-study.py', 'r') +print(f.read()) +f.close() + + for starttime in starttimes: for mesh_resolution, solver_tol in resolutions.items(): # initialise LDD simulation class