diff --git a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/Archive/TP-R-2-patch-mesh-study-new-gli-bak.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/Archive/TP-R-2-patch-mesh-study-new-gli-bak.py new file mode 100755 index 0000000000000000000000000000000000000000..23463c313e5d0a85d86406c0a734c1e3f7c389fd --- /dev/null +++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/Archive/TP-R-2-patch-mesh-study-new-gli-bak.py @@ -0,0 +1,512 @@ +#!/usr/bin/python3 +import dolfin as df +import mshr +import numpy as np +import sympy as sym +import typing as tp +import domainPatch as dp +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") +#import ufl as ufl + +# init sympy session +sym.init_printing() + +use_case = "TP-R-2-patch-realistic-new-nw-gli" +# solver_tol = 6E-7 +max_iter_num = 1000 +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, + } + +############ GRID ####################### +# mesh_resolution = 20 +timestep_size = 0.001 +number_of_timesteps = 800 +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] + +Lw = 0.025 #/timestep_size +Lnw=Lw + +lambda_w = 40 +lambda_nw = 40 + +include_gravity = False +debugflag = False +analyse_condition = False + +output_string = "./output/{}-{}_timesteps{}_P{}".format(datestr, use_case, number_of_timesteps, FEM_Lagrange_degree) + +# 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 +sub_domain0_vertices = [df.Point(-1.0, -1.0), + df.Point(1.0, -1.0), + df.Point(1.0, 1.0), + df.Point(-1.0, 1.0)] +# interface between subdomain1 and subdomain2 +interface12_vertices = [df.Point(-1.0, 0.0), + df.Point(1.0, 0.0) ] +# subdomain1. +sub_domain1_vertices = [interface12_vertices[0], + interface12_vertices[1], + sub_domain0_vertices[2], + sub_domain0_vertices[3]] + +# vertex coordinates of the outer boundaries. If it can not be specified as a +# polygon, use an entry per boundary polygon. This information is used for defining +# the Dirichlet boundary conditions. If a domain is completely internal, the +# dictionary entry should be 0: None +subdomain1_outer_boundary_verts = { + 0: [interface12_vertices[1], # + sub_domain0_vertices[2], + sub_domain0_vertices[3], # + interface12_vertices[0]] +} +# subdomain2 +sub_domain2_vertices = [sub_domain0_vertices[0], + sub_domain0_vertices[1], + interface12_vertices[1], + interface12_vertices[0] ] + +subdomain2_outer_boundary_verts = { + 0: [interface12_vertices[0], # + sub_domain0_vertices[0], + sub_domain0_vertices[1], + interface12_vertices[1]] +} + +# list of subdomains given by the boundary polygon vertices. +# Subdomains are given as a list of dolfin points forming +# a closed polygon, such that mshr.Polygon(subdomain_def_points[i]) can be used +# to create the subdomain. subdomain_def_points[0] contains the +# vertices of the global simulation domain and subdomain_def_points[i] contains the +# vertices of the subdomain i. +subdomain_def_points = [sub_domain0_vertices,# + sub_domain1_vertices,# + sub_domain2_vertices] +# in the below list, index 0 corresponds to the 12 interface which has index 1 +interface_def_points = [interface12_vertices] + +# if a subdomain has no outer boundary write None instead, i.e. +# i: None +# if i is the index of the inner subdomain. +outer_boundary_def_points = { + # subdomain number + 1 : subdomain1_outer_boundary_verts, + 2 : subdomain2_outer_boundary_verts +} + +# adjacent_subdomains[i] contains the indices of the subdomains sharing the +# interface i (i.e. given by interface_def_points[i]). +adjacent_subdomains = [[1,2]] +isRichards = { + 1: True, # + 2: False + } + + +viscosity = {# +# subdom_num : viscosity + 1 : {'wetting' :1}, + #'nonwetting': 1}, # + 2 : {'wetting' :1, + 'nonwetting': 1/50} +} + +porosity = {# +# subdom_num : porosity + 1 : 0.22,# + 2 : 0.0022 +} + +# Dict of the form: { subdom_num : density } +densities = { + 1: {'wetting': 997}, + 2: {'wetting': 997, + 'nonwetting': 1.225}, +} + +gravity_acceleration = 9.81 + +L = {# +# subdom_num : subdomain L for L-scheme + 1 : {'wetting' :Lw}, + # 'nonwetting': 0.25},# + 2 : {'wetting' :Lw, + 'nonwetting': Lnw} +} + + +lambda_param = {# +# subdom_num : lambda parameter for the L-scheme + 1 : {'wetting' :lambda_w, + 'nonwetting': lambda_nw},# + 2 : {'wetting' :lambda_w, + 'nonwetting': lambda_nw} +} + +## relative permeabilty functions on subdomain 1 +def rel_perm1w(s): + # relative permeabilty wetting on subdomain1 + return s**2 + +# def rel_perm1nw(s): +# # relative permeabilty nonwetting on subdomain1 +# return (1-s)**2 + +_rel_perm1w = ft.partial(rel_perm1w) +# _rel_perm1nw = ft.partial(rel_perm1nw) +subdomain1_rel_perm = { + 'wetting': _rel_perm1w,# + # 'nonwetting': _rel_perm1nw +} +## relative permeabilty functions on subdomain 2 +def rel_perm2w(s): + # relative permeabilty wetting on subdomain2 + return 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 + +_rel_perm2w = ft.partial(rel_perm2w) +_rel_perm2nw = ft.partial(rel_perm2nw) + +subdomain2_rel_perm = { + 'wetting': _rel_perm2w,# + 'nonwetting': _rel_perm2nw +} + +## dictionary of relative permeabilties on all domains. +relative_permeability = {# + 1: subdomain1_rel_perm, + 2: subdomain2_rel_perm +} + + +# definition of the derivatives of the relative permeabilities +# relative permeabilty functions on subdomain 1 +def rel_perm1w_prime(s): + # relative permeabilty on subdomain1 + return 2*s + +# def rel_perm1nw_prime(s): +# # relative permeabilty on subdomain1 +# return 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 + +def rel_perm2nw_prime(s): + # relative permeabilty on subdomain1 + return -3*(1-s)**2 + +_rel_perm1w_prime = ft.partial(rel_perm1w_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 +} + + +subdomain2_rel_perm_prime = { + 'wetting': _rel_perm2w_prime, + 'nonwetting': _rel_perm2nw_prime +} + +# dictionary of relative permeabilties on all domains. +ka_prime = { + 1: subdomain1_rel_perm_prime, + 2: subdomain2_rel_perm_prime, +} + + +# def saturation1(pc, subdomain_index): +# # inverse capillary pressure-saturation-relationship +# return df.conditional(pc > 0, 1/((1 + pc)**(1/(subdomain_index + 1))), 1) +# +# def saturation2(pc, n_index, alpha): +# # inverse capillary pressure-saturation-relationship +# return df.conditional(pc > 0, 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index)), 1) +# +# # S-pc-relation ship. We use the van Genuchten approach, i.e. pc = 1/alpha*(S^{-1/m} -1)^1/n, where +# # we set alpha = 0, assume m = 1-1/n (see Helmig) and assume that residual saturation is Sw +# def saturation1_sym(pc, subdomain_index): +# # inverse capillary pressure-saturation-relationship +# return 1/((1 + pc)**(1/(subdomain_index + 1))) +# +# +# def saturation2_sym(pc, n_index, alpha): +# # inverse capillary pressure-saturation-relationship +# #df.conditional(pc > 0, +# return 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index)) +# +# +# # derivative of S-pc relationship with respect to pc. This is needed for the +# # construction of a analytic solution. +# def saturation1_sym_prime(pc, subdomain_index): +# # inverse capillary pressure-saturation-relationship +# return -(1/(subdomain_index + 1))*(1 + pc)**((-subdomain_index - 2)/(subdomain_index + 1)) +# +# +# def saturation2_sym_prime(pc, n_index, alpha): +# # inverse capillary pressure-saturation-relationship +# return -(alpha*(n_index - 1)*(alpha*pc)**(n_index - 1)) / ( (1 + (alpha*pc)**n_index)**((2*n_index - 1)/n_index) ) +# +# # note that the conditional definition of S-pc in the nonsymbolic part will be +# # incorporated in the construction of the exact solution below. +# S_pc_sym = { +# 1: ft.partial(saturation1_sym, subdomain_index = 1), +# 2: ft.partial(saturation2_sym, n_index=3, alpha=0.001), +# } +# +# S_pc_sym_prime = { +# 1: ft.partial(saturation1_sym_prime, subdomain_index = 1), +# 2: ft.partial(saturation2_sym_prime, n_index=3, alpha=0.001), +# } +# +# sat_pressure_relationship = { +# 1: ft.partial(saturation1, subdomain_index = 1),#, +# 2: ft.partial(saturation2, n_index=3, alpha=0.001), +# } + +def saturation(pc, index): + # inverse capillary pressure-saturation-relationship + return df.conditional(pc > 0, 1/((1 + pc)**(1/(index + 1))), 1) + + +def saturation_sym(pc, index): + # inverse capillary pressure-saturation-relationship + return 1/((1 + pc)**(1/(index + 1))) + + +# derivative of S-pc relationship with respect to pc. This is needed for the +# construction of a analytic solution. +def saturation_sym_prime(pc, index): + # inverse capillary pressure-saturation-relationship + return -1/((index+1)*(1 + pc)**((index+2)/(index+1))) + + +# note that the conditional definition of S-pc in the nonsymbolic part will be +# incorporated in the construction of the exact solution below. +S_pc_sym = { + 1: ft.partial(saturation_sym, index=1), + 2: ft.partial(saturation_sym, index=2), + # 3: ft.partial(saturation_sym, index=2), + # 4: ft.partial(saturation_sym, index=1) +} + +S_pc_sym_prime = { + 1: ft.partial(saturation_sym_prime, index=1), + 2: ft.partial(saturation_sym_prime, index=2), + # 3: ft.partial(saturation_sym_prime, index=2), + # 4: ft.partial(saturation_sym_prime, index=1) +} + +sat_pressure_relationship = { + 1: ft.partial(saturation, index=1), + 2: ft.partial(saturation, index=2), + # 3: ft.partial(saturation, index=2), + # 4: ft.partial(saturation, index=1) +} + + +############################################# +# Manufacture source expressions with sympy # +############################################# +x, y = sym.symbols('x[0], x[1]') # needed by UFL +t = sym.symbols('t', positive=True) + +p_e_sym = { + 1: {'wetting': (-5.0 - (1.0 + t*t)*(1.0 + x*x + y*y))}, #*(1-x)**2*(1+x)**2*(1-y)**2}, + 2: {'wetting': (-5.0 - (1.0 + t*t)*(1.0 + x*x)), #*(1-x)**2*(1+x)**2*(1+y)**2, + 'nonwetting': (-1-t*(1.1+y + x**2))*y**3}, #*(1-x)**2*(1+x)**2*(1+y)**2}, +} #-y*y*(sym.sin(-2*t+2*x)*sym.sin(1/2*y-1.2*t)) - t*t*x*(0.5-y)*y*(1-x) + + +pc_e_sym = dict() +for subdomain, isR in isRichards.items(): + if isR: + pc_e_sym.update({subdomain: -p_e_sym[subdomain]['wetting'].copy()}) + else: + pc_e_sym.update({subdomain: p_e_sym[subdomain]['nonwetting'].copy() + - p_e_sym[subdomain]['wetting'].copy()}) + + +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, + ) +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() +# 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}} +for subdomain in isRichards.keys(): + # if subdomain has no outer boundary, outer_boundary_def_points[subdomain] is None + if outer_boundary_def_points[subdomain] is None: + dirichletBC.update({subdomain: None}) + else: + dirichletBC.update({subdomain: dict()}) + # set the dirichlet conditions to be the same code as exact solution on + # the subdomain. + for outer_boundary_ind in outer_boundary_def_points[subdomain].keys(): + dirichletBC[subdomain].update( + {outer_boundary_ind: exact_solution[subdomain]} + ) + + +# def saturation(pressure, subdomain_index): +# # inverse capillary pressure-saturation-relationship +# return df.conditional(pressure < 0, 1/((1 - pressure)**(1/(subdomain_index + 1))), 1) +# +# sa + +for starttime in starttimes: + 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-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/Archive/TP-R-2-patch-mesh-study-with-intrinsic-permeability.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/Archive/TP-R-2-patch-mesh-study-with-intrinsic-permeability.py new file mode 100755 index 0000000000000000000000000000000000000000..a253874797ad0e428c1badbc02ffef67d9341d5f --- /dev/null +++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/Archive/TP-R-2-patch-mesh-study-with-intrinsic-permeability.py @@ -0,0 +1,517 @@ +#!/usr/bin/python3 +import dolfin as df +import mshr +import numpy as np +import sympy as sym +import typing as tp +import domainPatch as dp +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") +#import ufl as ufl + +# init sympy session +sym.init_printing() + +use_case = "TPR-2-patch-realistic-with-intrinsic" +# solver_tol = 6E-7 +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, + } + +############ GRID ####################### +# mesh_resolution = 20 +timestep_size = 0.001 +number_of_timesteps = 800 +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] + +Lw = 0.025 #/timestep_size +Lnw=Lw + +lambda_w = 40 +lambda_nw = 40 + +include_gravity = False +debugflag = False +analyse_condition = True + +output_string = "./output/{}-{}_timesteps{}_P{}".format(datestr, use_case, number_of_timesteps, FEM_Lagrange_degree) + +# toggle what should be written to files +if mesh_study: + 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 + } +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 +sub_domain0_vertices = [df.Point(-1.0, -1.0), + df.Point(1.0, -1.0), + df.Point(1.0, 1.0), + df.Point(-1.0, 1.0)] +# interface between subdomain1 and subdomain2 +interface12_vertices = [df.Point(-1.0, 0.0), + df.Point(1.0, 0.0) ] +# subdomain1. +sub_domain1_vertices = [interface12_vertices[0], + interface12_vertices[1], + sub_domain0_vertices[2], + sub_domain0_vertices[3]] + +# vertex coordinates of the outer boundaries. If it can not be specified as a +# polygon, use an entry per boundary polygon. This information is used for defining +# the Dirichlet boundary conditions. If a domain is completely internal, the +# dictionary entry should be 0: None +subdomain1_outer_boundary_verts = { + 0: [interface12_vertices[1], # + sub_domain0_vertices[2], + sub_domain0_vertices[3], # + interface12_vertices[0]] +} +# subdomain2 +sub_domain2_vertices = [sub_domain0_vertices[0], + sub_domain0_vertices[1], + interface12_vertices[1], + interface12_vertices[0] ] + +subdomain2_outer_boundary_verts = { + 0: [interface12_vertices[0], # + sub_domain0_vertices[0], + sub_domain0_vertices[1], + interface12_vertices[1]] +} + +# list of subdomains given by the boundary polygon vertices. +# Subdomains are given as a list of dolfin points forming +# a closed polygon, such that mshr.Polygon(subdomain_def_points[i]) can be used +# to create the subdomain. subdomain_def_points[0] contains the +# vertices of the global simulation domain and subdomain_def_points[i] contains the +# vertices of the subdomain i. +subdomain_def_points = [sub_domain0_vertices,# + sub_domain1_vertices,# + sub_domain2_vertices] +# in the below list, index 0 corresponds to the 12 interface which has index 1 +interface_def_points = [interface12_vertices] + +# if a subdomain has no outer boundary write None instead, i.e. +# i: None +# if i is the index of the inner subdomain. +outer_boundary_def_points = { + # subdomain number + 1 : subdomain1_outer_boundary_verts, + 2 : subdomain2_outer_boundary_verts +} + +# adjacent_subdomains[i] contains the indices of the subdomains sharing the +# interface i (i.e. given by interface_def_points[i]). +adjacent_subdomains = [[1,2]] +isRichards = { + 1: True, # + 2: False + } + + +viscosity = {# +# subdom_num : viscosity + 1 : {'wetting' :1}, + #'nonwetting': 1}, # + 2 : {'wetting' :1, + 'nonwetting': 1/50} +} + +porosity = {# +# subdom_num : porosity + 1 : 0.22,# + 2 : 0.0022 +} + +# Dict of the form: { subdom_num : density } +densities = { + 1: {'wetting': 997}, + 2: {'wetting': 997, + 'nonwetting': 1.225}, +} + +intrinsic_permeability = { + 1: 10e-11, #laut wikipedia in der range of sand + 2: 10e-14, # layered clay achtung: auf wikipedia in cm^2 angegeb. hier m^2 +} + +gravity_acceleration = 9.81 + +L = {# +# subdom_num : subdomain L for L-scheme + 1 : {'wetting' :Lw}, + # 'nonwetting': 0.25},# + 2 : {'wetting' :Lw, + 'nonwetting': Lnw} +} + + +lambda_param = {# +# subdom_num : lambda parameter for the L-scheme + 1 : {'wetting' :lambda_w, + 'nonwetting': lambda_nw},# + 2 : {'wetting' :lambda_w, + 'nonwetting': lambda_nw} +} + +## relative permeabilty functions on subdomain 1 +def rel_perm1w(s): + # relative permeabilty wetting on subdomain1 + return intrinsic_permeability[1]*s**2 + +# def rel_perm1nw(s): +# # relative permeabilty nonwetting on subdomain1 +# return (1-s)**2 + +_rel_perm1w = ft.partial(rel_perm1w) +# _rel_perm1nw = ft.partial(rel_perm1nw) +subdomain1_rel_perm = { + 'wetting': _rel_perm1w,# + # 'nonwetting': _rel_perm1nw +} +## relative permeabilty functions on subdomain 2 +def rel_perm2w(s): + # relative permeabilty wetting on subdomain2 + 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 intrinsic_permeability[1]*(1-s)**3 + +_rel_perm2w = ft.partial(rel_perm2w) +_rel_perm2nw = ft.partial(rel_perm2nw) + +subdomain2_rel_perm = { + 'wetting': _rel_perm2w,# + 'nonwetting': _rel_perm2nw +} + +## dictionary of relative permeabilties on all domains. +relative_permeability = {# + 1: subdomain1_rel_perm, + 2: subdomain2_rel_perm +} + + +# definition of the derivatives of the relative permeabilities +# relative permeabilty functions on subdomain 1 +def rel_perm1w_prime(s): + # relative permeabilty on subdomain1 + return 2*s + +# def rel_perm1nw_prime(s): +# # relative permeabilty on subdomain1 +# return 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 + +def rel_perm2nw_prime(s): + # relative permeabilty on subdomain1 + return -3*(1-s)**2 + +_rel_perm1w_prime = ft.partial(rel_perm1w_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 +} + + +subdomain2_rel_perm_prime = { + 'wetting': _rel_perm2w_prime, + 'nonwetting': _rel_perm2nw_prime +} + +# dictionary of relative permeabilties on all domains. +ka_prime = { + 1: subdomain1_rel_perm_prime, + 2: subdomain2_rel_perm_prime, +} + + +# def saturation1(pc, subdomain_index): +# # inverse capillary pressure-saturation-relationship +# return df.conditional(pc > 0, 1/((1 + pc)**(1/(subdomain_index + 1))), 1) +# +# def saturation2(pc, n_index, alpha): +# # inverse capillary pressure-saturation-relationship +# return df.conditional(pc > 0, 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index)), 1) +# +# # S-pc-relation ship. We use the van Genuchten approach, i.e. pc = 1/alpha*(S^{-1/m} -1)^1/n, where +# # we set alpha = 0, assume m = 1-1/n (see Helmig) and assume that residual saturation is Sw +# def saturation1_sym(pc, subdomain_index): +# # inverse capillary pressure-saturation-relationship +# return 1/((1 + pc)**(1/(subdomain_index + 1))) +# +# +# def saturation2_sym(pc, n_index, alpha): +# # inverse capillary pressure-saturation-relationship +# #df.conditional(pc > 0, +# return 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index)) +# +# +# # derivative of S-pc relationship with respect to pc. This is needed for the +# # construction of a analytic solution. +# def saturation1_sym_prime(pc, subdomain_index): +# # inverse capillary pressure-saturation-relationship +# return -(1/(subdomain_index + 1))*(1 + pc)**((-subdomain_index - 2)/(subdomain_index + 1)) +# +# +# def saturation2_sym_prime(pc, n_index, alpha): +# # inverse capillary pressure-saturation-relationship +# return -(alpha*(n_index - 1)*(alpha*pc)**(n_index - 1)) / ( (1 + (alpha*pc)**n_index)**((2*n_index - 1)/n_index) ) +# +# # note that the conditional definition of S-pc in the nonsymbolic part will be +# # incorporated in the construction of the exact solution below. +# S_pc_sym = { +# 1: ft.partial(saturation1_sym, subdomain_index = 1), +# 2: ft.partial(saturation2_sym, n_index=3, alpha=0.001), +# } +# +# S_pc_sym_prime = { +# 1: ft.partial(saturation1_sym_prime, subdomain_index = 1), +# 2: ft.partial(saturation2_sym_prime, n_index=3, alpha=0.001), +# } +# +# sat_pressure_relationship = { +# 1: ft.partial(saturation1, subdomain_index = 1),#, +# 2: ft.partial(saturation2, n_index=3, alpha=0.001), +# } + +def saturation(pc, index): + # inverse capillary pressure-saturation-relationship + return df.conditional(pc > 0, 1/((1 + pc)**(1/(index + 1))), 1) + + +def saturation_sym(pc, index): + # inverse capillary pressure-saturation-relationship + return 1/((1 + pc)**(1/(index + 1))) + + +# derivative of S-pc relationship with respect to pc. This is needed for the +# construction of a analytic solution. +def saturation_sym_prime(pc, index): + # inverse capillary pressure-saturation-relationship + return -1/((index+1)*(1 + pc)**((index+2)/(index+1))) + + +# note that the conditional definition of S-pc in the nonsymbolic part will be +# incorporated in the construction of the exact solution below. +S_pc_sym = { + 1: ft.partial(saturation_sym, index=1), + 2: ft.partial(saturation_sym, index=2), + # 3: ft.partial(saturation_sym, index=2), + # 4: ft.partial(saturation_sym, index=1) +} + +S_pc_sym_prime = { + 1: ft.partial(saturation_sym_prime, index=1), + 2: ft.partial(saturation_sym_prime, index=2), + # 3: ft.partial(saturation_sym_prime, index=2), + # 4: ft.partial(saturation_sym_prime, index=1) +} + +sat_pressure_relationship = { + 1: ft.partial(saturation, index=1), + 2: ft.partial(saturation, index=2), + # 3: ft.partial(saturation, index=2), + # 4: ft.partial(saturation, index=1) +} + + +############################################# +# Manufacture source expressions with sympy # +############################################# +x, y = sym.symbols('x[0], x[1]') # needed by UFL +t = sym.symbols('t', positive=True) + +p_e_sym = { + 1: {'wetting': (-6.0 - (1.0 + t*t)*(1.0 + x*x + y*y))}, #*(1-x)**2*(1+x)**2*(1-y)**2}, + 2: {'wetting': (-6.0 - (1.0 + t*t)*(1.0 + x*x)), #*(1-x)**2*(1+x)**2*(1+y)**2, + 'nonwetting': (-1-t*(1.1+y + x**2))*y**2}, #*(1-x)**2*(1+x)**2*(1+y)**2}, +} #-y*y*(sym.sin(-2*t+2*x)*sym.sin(1/2*y-1.2*t)) - t*t*x*(0.5-y)*y*(1-x) + + +pc_e_sym = dict() +for subdomain, isR in isRichards.items(): + if isR: + pc_e_sym.update({subdomain: -p_e_sym[subdomain]['wetting'].copy()}) + else: + pc_e_sym.update({subdomain: p_e_sym[subdomain]['nonwetting'].copy() + - p_e_sym[subdomain]['wetting'].copy()}) + + +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, + ) +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() +# 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}} +for subdomain in isRichards.keys(): + # if subdomain has no outer boundary, outer_boundary_def_points[subdomain] is None + if outer_boundary_def_points[subdomain] is None: + dirichletBC.update({subdomain: None}) + else: + dirichletBC.update({subdomain: dict()}) + # set the dirichlet conditions to be the same code as exact solution on + # the subdomain. + for outer_boundary_ind in outer_boundary_def_points[subdomain].keys(): + dirichletBC[subdomain].update( + {outer_boundary_ind: exact_solution[subdomain]} + ) + + +# def saturation(pressure, subdomain_index): +# # inverse capillary pressure-saturation-relationship +# return df.conditional(pressure < 0, 1/((1 - pressure)**(1/(subdomain_index + 1))), 1) +# +# sa + +for starttime in starttimes: + 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-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-realistic-no-gravity-but-varying-intrinsic-perm-mesh-study.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/Archive/TP-R-2-patch-realistic-no-gravity-but-varying-intrinsic-perm-mesh-study.py similarity index 100% rename from Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-realistic-no-gravity-but-varying-intrinsic-perm-mesh-study.py rename to Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/Archive/TP-R-2-patch-realistic-no-gravity-but-varying-intrinsic-perm-mesh-study.py diff --git a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-realistic-no-gravity-but-varying-intrinsic-perm.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-realistic-no-gravity-but-varying-intrinsic-perm.py new file mode 100755 index 0000000000000000000000000000000000000000..280a0fb59b887d6e0efe14922cb8f237d39dbe44 --- /dev/null +++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-realistic-no-gravity-but-varying-intrinsic-perm.py @@ -0,0 +1,343 @@ +#!/usr/bin/python3 +"""TPR 2 patch soil simulation. + +This program sets up an LDD simulation +""" +import dolfin as df +import sympy as sym +import functions as fts +import LDDsimulation as ldd +import helpers as hlp +import datetime +import os +import multiprocessing as mp +import domainSubstructuring as dss + +# init sympy session +sym.init_printing() + +# PREREQUISITS ############################################################### +# check if output directory "./output" exists. This will be used in +# the generation of the output string. +if not os.path.exists('./output'): + os.mkdir('./output') + print("Directory ", './output', " created ") +else: + print("Directory ", './output', " already exists. Will use as output \ + directory") + +date = datetime.datetime.now() +datestr = date.strftime("%Y-%m-%d") + +# Name of the usecase that will be printed during simulation. +use_case = "TP-R-2-patch-realistic-no-gravity-but-varying-intrinsic-perm-flatL" +# The name of this very file. Needed for creating log output. +thisfile = "TP-R-2-patch-realistic-no-gravity-but-varying-intrinsic-perm.py" + +# GENERAL SOLVER CONFIG ###################################################### +# maximal iteration per timestep +max_iter_num = 1000 +FEM_Lagrange_degree = 1 + +# GRID AND MESH STUDY SPECIFICATIONS ######################################### +mesh_study = True +resolutions = { + 1: 1e-5, + 2: 1e-5, + 4: 1e-5, + 8: 1e-5, + 16: 3e-6, + 32: 1e-6, + 64: 1e-6, + 128: 1e-6, + # 256: 1e-6, + } + +# starttimes gives a list of starttimes to run the simulation from. +# The list is looped over and a simulation is run with t_0 as initial time +# for each element t_0 in starttimes. +# starttimes = {0: 0.0, 1: 0.3, 2: 0.6, 3: 0.9} +starttimes = {0: 0.0} +timestep_size = 0.01 +number_of_timesteps = 1000 + +# LDD scheme parameters ###################################################### +Lw1 = 0.007 +Lnw1 = 0.007 + +Lw2 = 0.0005 +Lnw2 = 0.0005 + +lambda_w = 0.01 +lambda_nw = 0.01 + +include_gravity = False +debugflag = False +analyse_condition = False + +# I/O CONFIG ################################################################# +# when number_of_timesteps is high, it might take a long time to write all +# timesteps to disk. Therefore, you can choose to only write data of every +# plot_timestep_every timestep to disk. +plot_timestep_every = 3 +# Decide how many timesteps you want analysed. Analysed means, that +# subsequent errors of the L-iteration within the timestep are written out. +number_of_timesteps_to_analyse = 10 + +# fine grained control over data to be written to disk in the mesh study case +# as well as for a regular simuation for a fixed grid. +if mesh_study: + write_to_file = { + # output the relative errornorm (integration in space) w.r.t. an exact + # solution for each timestep into a csv file. + 'space_errornorms': True, + # save the mesh and marker functions to disk + 'meshes_and_markers': True, + # save xdmf/h5 data for each LDD iteration for timesteps determined by + # number_of_timesteps_to_analyse. I/O intensive! + 'L_iterations_per_timestep': False, + # save solution to xdmf/h5. + 'solutions': True, + # save absolute differences w.r.t an exact solution to xdmf/h5 file + # to monitor where on the domains errors happen + 'absolute_differences': True, + # analyise condition numbers for timesteps determined by + # number_of_timesteps_to_analyse and save them over time to csv. + 'condition_numbers': analyse_condition, + # output subsequent iteration errors measured in L^2 to csv for + # timesteps determined by number_of_timesteps_to_analyse. + # Usefull to monitor convergence of the acutal LDD solver. + 'subsequent_errors': True + } +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 + } + +# OUTPUT FILE STRING ######################################################### +output_string = "./output/{}-{}_timesteps{}_P{}".format( + datestr, use_case, number_of_timesteps, FEM_Lagrange_degree + ) + +# DOMAIN AND INTERFACE ####################################################### +substructuring = dss.twoSoilLayers() +interface_def_points = substructuring.interface_def_points +adjacent_subdomains = substructuring.adjacent_subdomains +subdomain_def_points = substructuring.subdomain_def_points +outer_boundary_def_points = substructuring.outer_boundary_def_points + +# MODEL CONFIGURATION ######################################################### +isRichards = { + 1: True, + 2: False + } + + +viscosity = { + # subdom_num : viscosity + 1: {'wetting': 1, + 'nonwetting': 1/50}, + 2: {'wetting': 1, + 'nonwetting': 1/50} +} + +porosity = { + # subdom_num : porosity + 1: 0.22, + 2: 0.022 +} + +# Dict of the form: { subdom_num : density } +densities = { + 1: {'wetting': 997, + 'nonwetting': 1.225}, + 2: {'wetting': 997, + 'nonwetting': 1.225} +} + +gravity_acceleration = 9.81 + +L = { + # subdom_num : subdomain L for L-scheme + 1: {'wetting': Lw1, + 'nonwetting': Lnw1}, + 2: {'wetting': Lw2, + 'nonwetting': Lnw2} +} + + +lambda_param = { + # interface_num : lambda parameter for the L-scheme + 0: {'wetting': lambda_w, + 'nonwetting': lambda_nw}, +} + +intrinsic_permeability = { + 1: 0.01, + 2: 0.0001, +} + +# relative permeabilties +rel_perm_definition = { + 1: {"wetting": "Spow2", + "nonwetting": "oneMinusSpow2"}, + 2: {"wetting": "Spow3", + "nonwetting": "oneMinusSpow3"}, +} + +rel_perm_dict = fts.generate_relative_permeability_dicts(rel_perm_definition) +relative_permeability = rel_perm_dict["ka"] +ka_prime = rel_perm_dict["ka_prime"] + +# S-pc relation +Spc_on_subdomains = { + 1: {"testSpc": {"index": 1}}, + 2: {"testSpc": {"index": 2}}, +} + +Spc = fts.generate_Spc_dicts(Spc_on_subdomains) +S_pc_sym = Spc["symbolic"] +S_pc_sym_prime = Spc["prime_symbolic"] +sat_pressure_relationship = Spc["dolfin"] + +############################################################################### +# Manufacture source expressions with sympy # +############################################################################### +x, y = sym.symbols('x[0], x[1]') # needed by UFL +t = sym.symbols('t', positive=True) + +p_e_sym = { + 1: {'wetting': (-7.0 - (1.0 + t*t)*(1.0 + x*x + y*y))}, #*(1-x)**2*(1+x)**2*(1-y)**2}, + 2: {'wetting': (-7.0 - (1.0 + t*t)*(1.0 + x*x)), #*(1-x)**2*(1+x)**2*(1+y)**2, + 'nonwetting': (-2-t*(1.1+y + x**2))*y**2}, #*(1-x)**2*(1+x)**2*(1+y)**2}, +} #-y*y*(sym.sin(-2*t+2*x)*sym.sin(1/2*y-1.2*t)) - t*t*x*(0.5-y)*y*(1-x) + +pc_e_sym = hlp.generate_exact_symbolic_pc( + isRichards=isRichards, + symbolic_pressure=p_e_sym + ) + +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, + intrinsic_permeability=intrinsic_permeability, + relative_permeability=relative_permeability, + relative_permeability_prime=ka_prime, + densities=densities, + gravity_acceleration=gravity_acceleration, + include_gravity=include_gravity, + ) +source_expression = exact_solution_example['source'] +exact_solution = exact_solution_example['exact_solution'] +initial_condition = exact_solution_example['initial_condition'] + +# BOUNDARY CONDITIONS ######################################################### +# Dictionary of dirichlet boundary conditions. If an exact solution case is +# used, use the hlp.generate_exact_DirichletBC() method to generate the +# Dirichlet Boundary conditions from the exact solution. +dirichletBC = hlp.generate_exact_DirichletBC( + isRichards=isRichards, + outer_boundary_def_points=outer_boundary_def_points, + exact_solution=exact_solution + ) +# If no exact solution is provided you need to provide a dictionary of boundary +# conditions. See the definiton of hlp.generate_exact_DirichletBC() to see +# the structure. + +# LOG FILE OUTPUT ############################################################# +# read this file and print it to std out. This way the simulation can produce a +# log file with ./TP-R-layered_soil.py | tee simulation.log +f = open(thisfile, 'r') +print(f.read()) +f.close() + +# MAIN ######################################################################## +if __name__ == '__main__': + # dictionary of simualation parameters to pass to the run function. + # mesh_resolution and starttime are excluded, as they get passed explicitly + # to achieve parallelisation in these parameters in these parameters for + # mesh studies etc. + simulation_parameter = { + "tol": 1E-14, + "debugflag": debugflag, + "max_iter_num": max_iter_num, + "FEM_Lagrange_degree": FEM_Lagrange_degree, + "mesh_study": mesh_study, + "use_case": use_case, + "output_string": 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, + "intrinsic_permeability": intrinsic_permeability, + "sat_pressure_relationship": 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, + "source_expression": source_expression, + "initial_condition": initial_condition, + "dirichletBC": dirichletBC, + "exact_solution": exact_solution, + "densities": densities, + "include_gravity": include_gravity, + "gravity_acceleration": gravity_acceleration, + "write_to_file": write_to_file, + "analyse_condition": analyse_condition + } + for number_shift, starttime in starttimes.items(): + simulation_parameter.update( + {"starttime_timestep_number_shift": number_shift} + ) + for mesh_resolution, solver_tol in resolutions.items(): + simulation_parameter.update({"solver_tol": solver_tol}) + hlp.info(simulation_parameter["use_case"]) + processQueue = mp.Queue() + LDDsim = mp.Process( + target=hlp.run_simulation, + args=( + simulation_parameter, + processQueue, + starttime, + mesh_resolution + ) + ) + LDDsim.start() + # LDDsim.join() + # hlp.run_simulation( + # mesh_resolution=mesh_resolution, + # starttime=starttime, + # parameter=simulation_parameter + # ) + + # LDDsim.join() + if mesh_study: + simulation_output_dir = processQueue.get() + hlp.merge_spacetime_errornorms(isRichards=isRichards, + resolutions=resolutions, + output_dir=simulation_output_dir) diff --git a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-realistic-same-intrinsic-perm.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-realistic-same-intrinsic-perm.py index 6943cfe34cbb32e8ce1f56c00255fdf8da8903c5..3cb92069c56ce4412cbb07432db11562853895a1 100755 --- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-realistic-same-intrinsic-perm.py +++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-realistic-same-intrinsic-perm.py @@ -30,13 +30,13 @@ date = datetime.datetime.now() datestr = date.strftime("%Y-%m-%d") # Name of the usecase that will be printed during simulation. -use_case = "TP-R-2-patch-realistic-same-intrinsic-perm" +use_case = "TP-R-2-patch-realistic-same-intrinsic-better-L" # The name of this very file. Needed for creating log output. thisfile = "TP-R-2-patch-realistic-same-intrinsic-perm.py" # GENERAL SOLVER CONFIG ###################################################### # maximal iteration per timestep -max_iter_num = 250 +max_iter_num = 1000 FEM_Lagrange_degree = 1 # GRID AND MESH STUDY SPECIFICATIONS ######################################### @@ -46,29 +46,30 @@ resolutions = { 2: 1e-5, 4: 1e-5, 8: 1e-5, - 16: 1e-5, - 32: 1e-5, - 64: 5e-6, - 128: 5e-6, + 16: 3e-6, + 32: 1e-6, + 64: 1e-6, + 128: 1e-6, # 256: 1e-6, } # starttimes gives a list of starttimes to run the simulation from. # The list is looped over and a simulation is run with t_0 as initial time # for each element t_0 in starttimes. +# starttimes = {0: 0.0, 1: 0.3, 2: 0.6, 3: 0.9} starttimes = {0: 0.0} timestep_size = 0.001 -number_of_timesteps = 800 +number_of_timesteps = 1000 # LDD scheme parameters ###################################################### -Lw1 = 0.25 -Lnw1 = 0.25 +Lw1 = 0.007 +Lnw1 = 0.007 -Lw2 = 0.5 -Lnw2 = 0.25 +Lw2 = 0.005 +Lnw2 = 0.005 -lambda_w = 40 -lambda_nw = 40 +lambda_w = 0.75 +lambda_nw = 0.75 include_gravity = False debugflag = False @@ -78,10 +79,10 @@ analyse_condition = False # when number_of_timesteps is high, it might take a long time to write all # timesteps to disk. Therefore, you can choose to only write data of every # plot_timestep_every timestep to disk. -plot_timestep_every = 4 +plot_timestep_every = 3 # Decide how many timesteps you want analysed. Analysed means, that # subsequent errors of the L-iteration within the timestep are written out. -number_of_timesteps_to_analyse = 5 +number_of_timesteps_to_analyse = 10 # fine grained control over data to be written to disk in the mesh study case # as well as for a regular simuation for a fixed grid. @@ -140,9 +141,9 @@ isRichards = { viscosity = { # subdom_num : viscosity - 1: {'wetting': 1, + 1: {'wetting': 1.0, 'nonwetting': 1/50}, - 2: {'wetting': 1, + 2: {'wetting': 1.0, 'nonwetting': 1/50} } @@ -182,7 +183,7 @@ intrinsic_permeability = { 2: 0.01, } -# RELATIVE PEMRMEABILITIES +# relative permeabilties rel_perm_definition = { 1: {"wetting": "Spow2", "nonwetting": "oneMinusSpow2"},