diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py index 340b08c1309c69c82374247f1be9d0916b8cd4dc..ae785b2fcd89d4f9beebf7fee11614f0a9c4a0bc 100644 --- a/LDDsimulation/LDDsimulation.py +++ b/LDDsimulation/LDDsimulation.py @@ -265,6 +265,7 @@ class LDDsimulation(object): raise RuntimeError("include_gravity = True but no densities were given.") self._parameters_set = True + def initialise(self) -> None: """ initialise LDD simulation diff --git a/LDDsimulation/domainPatch.py b/LDDsimulation/domainPatch.py index 0b11785b48f362ca4e85c24eb0f41666b7aec23e..52476d16608b694b12a1f4d44f23ae85174e62dc 100644 --- a/LDDsimulation/domainPatch.py +++ b/LDDsimulation/domainPatch.py @@ -621,7 +621,7 @@ class DomainPatch(df.SubDomain): # print(f"np.allclose({p_coord}, {gli_dof_coord}, rtol=1e-10, atol=1e-14): {b_close_a}") if a_close_b and b_close_a: p_previous_dof = p_prev_dof - print("\n") + if p_previous_dof is None: raise RuntimeError(f"didn't find p_previous_dof corresponding to dict key ({gli_dof_coord})", "something is wrong") diff --git a/RR-2-patch-test-case/RR-2-patch-test.py b/RR-2-patch-test-case/RR-2-patch-test.py index 343a5b0c4edfea1df528d1d9de1208da5797e364..65e7f5340494d7845878e506ca461f7a9df34915 100755 --- a/RR-2-patch-test-case/RR-2-patch-test.py +++ b/RR-2-patch-test-case/RR-2-patch-test.py @@ -15,8 +15,8 @@ solver_tol = 1E-6 ############ GRID #######################ü mesh_resolution = 40 -timestep_size = 0.01 -number_of_timesteps = 15 +timestep_size = 0.001 +number_of_timesteps = 1500 # 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 @@ -25,11 +25,11 @@ starttime = 0 Lw = 0.25 #/timestep_size Lnw=Lw -l_param_w = 4 +lambda_w = 4 -include_gravity = False +include_gravity = True debugflag = False -analyse_condition = False +analyse_condition = True output_string = "./output/new_gravity_term-number_of_timesteps{}_".format(number_of_timesteps) @@ -138,8 +138,8 @@ L = {# lambda_param = {# # subdom_num : lambda parameter for the L-scheme - 1 : {'wetting': l_param_w},# - 2 : {'wetting': l_param_w} + 1 : {'wetting': lambda_w},# + 2 : {'wetting': lambda_w} } ## relative permeabilty functions on subdomain 1 @@ -226,7 +226,7 @@ write_to_file = { } # initialise LDD simulation class -simulation = ldd.LDDsimulation(tol = 1E-14, debug = False, LDDsolver_tol=1E-9) +simulation = ldd.LDDsimulation(tol = 1E-14, debug = debugflag, LDDsolver_tol=solver_tol) simulation.set_parameters(use_case=use_case, output_dir = output_string,# subdomain_def_points = subdomain_def_points,# @@ -255,7 +255,7 @@ simulation.set_parameters(use_case=use_case, ) simulation.initialise() -print(simulation.__dict__) +#print(simulation.__dict__) simulation.run(analyse_condition=analyse_condition) # simulation.LDDsolver(time = 0, debug = True, analyse_timestep = True) # df.info(parameters, True) diff --git a/RR-multi-patch-with-inner-patch/RR-multi-patch-with-inner-patch.py b/RR-multi-patch-with-inner-patch/RR-multi-patch-with-inner-patch.py index d83d619c36bfbe41f07718d182690f46e127f56c..df7dd19e9af0a15f802a322ada790241bd9a308d 100755 --- a/RR-multi-patch-with-inner-patch/RR-multi-patch-with-inner-patch.py +++ b/RR-multi-patch-with-inner-patch/RR-multi-patch-with-inner-patch.py @@ -7,28 +7,34 @@ import sympy as sym # import domainPatch as dp import LDDsimulation as ldd import functools as ft +import helpers as hlp # import ufl as ufl # init sympy session sym.init_printing() -# ----------------------------------------------------------------------------# -# ------------------- MESH ---------------------------------------------------# -# ----------------------------------------------------------------------------# -mesh_resolution = 51 -# ----------------------------------------:-------------------------------------# -# ------------------- TIME ---------------------------------------------------# -# ----------------------------------------------------------------------------# +use_case = "RR-multi-patch-with-inner-patch-cuttoff" +solver_tol = 1E-6 + +############ GRID #######################ü +mesh_resolution = 20 timestep_size = 0.01 number_of_timesteps = 150 -# decide how many timesteps you want analysed. Analysed means, that we write -# out subsequent errors of the L-iteration within the timestep. +# 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 starttime = 0 -# Lw = 0.25/timestep_size -Lw = 0.25 +Lw = 0.25 #/timestep_size +Lnw=Lw + +lambda_w = 4 +include_gravity = True +debugflag = False +analyse_condition = False + +output_string = "./output/new_gravity_cutoff_term-number_of_timesteps{}_".format(number_of_timesteps) # ----------------------------------------------------------------------------# # ------------------- Domain and Interface -----------------------------------# @@ -205,15 +211,15 @@ L = { 5: {'wetting': Lw} } -lamdal_w = 32 +lambda_w = 32 # subdom_num : lambda parameter for the L-scheme lambda_param = { - 1: {'wetting': lamdal_w}, - 2: {'wetting': lamdal_w}, - 3: {'wetting': lamdal_w}, - 4: {'wetting': lamdal_w}, - 5: {'wetting': lamdal_w} + 1: {'wetting': lambda_w}, + 2: {'wetting': lambda_w}, + 3: {'wetting': lambda_w}, + 4: {'wetting': lambda_w}, + 5: {'wetting': lambda_w} } @@ -460,84 +466,37 @@ p_e_sym = { } -pc_e_sym = { - 1: -1*p_e_sym[1]['wetting'], - 2: -1*p_e_sym[2]['wetting'], - 3: -1*p_e_sym[3]['wetting'], - 4: -1*p_e_sym[4]['wetting'], - 5: -1*p_e_sym[5]['wetting'] -} - -# turn above symbolic code into exact solution for dolphin and -# construct the rhs that matches the above exact solution. -dtS = dict() -div_flux = dict() -source_expression = dict() -exact_solution = dict() -initial_condition = dict() +pc_e_sym = dict() for subdomain, isR in isRichards.items(): - dtS.update({subdomain: dict()}) - div_flux.update({subdomain: dict()}) - source_expression.update({subdomain: dict()}) - exact_solution.update({subdomain: dict()}) - initial_condition.update({subdomain: dict()}) if isR: - subdomain_has_phases = ["wetting"] + pc_e_sym.update({subdomain: -p_e_sym[subdomain]['wetting']}) else: - subdomain_has_phases = ["wetting", "nonwetting"] - - # conditional for S_pc_prime - pc = pc_e_sym[subdomain] - dtpc = sym.diff(pc, t, 1) - dxpc = sym.diff(pc, x, 1) - dypc = sym.diff(pc, y, 1) - S = sym.Piecewise((S_pc_sym[subdomain](pc), pc > 0), (1, True)) - dS = sym.Piecewise((S_pc_sym_prime[subdomain](pc), pc > 0), (0, True)) - for phase in subdomain_has_phases: - # Turn above symbolic expression for exact solution into c code - exact_solution[subdomain].update( - {phase: sym.printing.ccode(p_e_sym[subdomain][phase])} - ) - # save the c code for initial conditions - initial_condition[subdomain].update( - {phase: sym.printing.ccode(p_e_sym[subdomain][phase].subs(t, 0))} - ) - if phase == "nonwetting": - dtS[subdomain].update( - {phase: -porosity[subdomain]*dS*dtpc} - ) - else: - dtS[subdomain].update( - {phase: porosity[subdomain]*dS*dtpc} - ) - pa = p_e_sym[subdomain][phase] - dxpa = sym.diff(pa, x, 1) - dxdxpa = sym.diff(pa, x, 2) - dypa = sym.diff(pa, y, 1) - dydypa = sym.diff(pa, y, 2) - mu = viscosity[subdomain][phase] - ka = relative_permeability[subdomain][phase] - dka = ka_prime[subdomain][phase] - rho = densities[subdomain][phase] - g = gravity_acceleration - - if phase == "nonwetting": - # x part of div(flux) for nonwetting - dxdxflux = -1/mu*dka(1-S)*dS*dxpc*dxpa + 1/mu*dxdxpa*ka(1-S) - # y part of div(flux) for nonwetting - dydyflux = -1/mu*dka(1-S)*dS*dypc*(dypa - rho*g) \ - + 1/mu*dydypa*ka(1-S) - else: - # x part of div(flux) for wetting - dxdxflux = 1/mu*dka(S)*dS*dxpc*dxpa + 1/mu*dxdxpa*ka(S) - # y part of div(flux) for wetting - dydyflux = 1/mu*dka(S)*dS*dypc*(dypa - rho*g) + 1/mu*dydypa*ka(S) - div_flux[subdomain].update({phase: dxdxflux + dydyflux}) - contructed_rhs = dtS[subdomain][phase] - div_flux[subdomain][phase] - source_expression[subdomain].update( - {phase: sym.printing.ccode(contructed_rhs)} - ) - # print(f"source_expression[{subdomain}][{phase}] =", source_expression[subdomain][phase]) + pc_e_sym.update({subdomain: p_e_sym[subdomain]['nonwetting'] + - p_e_sym[subdomain]['wetting']}) + +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() @@ -568,87 +527,15 @@ for subdomain in isRichards.keys(): ) -# # construction of the rhs that matches the above exact solution. -# dtS = dict() -# div_flux = dict() -# source_expression = dict() -# for subdomain, isR in isRichards.items(): -# dtS.update({subdomain: dict()}) -# div_flux.update({subdomain: dict()}) -# source_expression.update({subdomain: dict()}) -# if isR: -# subdomain_has_phases = ["wetting"] -# else: -# subdomain_has_phases = ["wetting", "nonwetting"] -# -# # conditional for S_pc_prime -# pc = pc_e_sym[subdomain] -# dtpc = sym.diff(pc, t, 1) -# dxpc = sym.diff(pc, x, 1) -# dypc = sym.diff(pc, y, 1) -# S = sym.Piecewise((S_pc_sym[subdomain](pc), pc > 0), (1, True)) -# dS = sym.Piecewise((S_pc_sym_prime[subdomain](pc), pc > 0), (0, True)) -# for phase in subdomain_has_phases: -# if phase == "nonwetting": -# dS = -dS -# dtS[subdomain].update( -# {phase: porosity[subdomain]*dS*dtpc} -# ) -# pa = p_e_sym[subdomain][phase] -# dxpa = sym.diff(pa, x, 1) -# dxdxpa = sym.diff(pa, x, 2) -# dypa = sym.diff(pa, y, 1) -# dydypa = sym.diff(pa, y, 2) -# mu = viscosity[subdomain][phase] -# ka = relative_permeability[subdomain][phase] -# dka = ka_prime[subdomain][phase] -# rho = densities[subdomain][phase] -# g = gravity_acceleration -# -# if phase == "nonwetting": -# dxdxflux = 1/mu*dka(1-S)*dS*dxpc*dxpa - 1/mu*dxdxpa*ka(1-S) -# dydyflux = 1/mu*dka(1-S)*dS*dypc*(dypa - rho*g) \ -# - 1/mu*dydypa*ka(1-S) -# else: -# dxdxflux = -1/mu*dka(S)*dS*dxpc*dxpa - 1/mu*dxdxpa*ka(S) -# dydyflux = -1/mu*dka(S)*dS*dypc*(dypa - rho*g) - 1/mu*dydypa*ka(S) -# div_flux[subdomain].update({phase: dxdxflux + dydyflux}) -# contructed_rhs = dtS[subdomain][phase] + div_flux[subdomain][phase] -# source_expression[subdomain].update( -# {phase: sym.printing.ccode(contructed_rhs)} -# ) -# print(f"source_expression[{subdomain}][{phase}] =", -# source_expression[subdomain][phase]) - -# 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}} -# dirichletBC = { -# 1: {0: {'wetting': exact_solution[1]['wetting']}}, -# 2: {0: {'wetting': exact_solution[2]['wetting']}}, -# 3: {0: {'wetting': exact_solution[3]['wetting']}}, -# 4: {0: {'wetting': exact_solution[4]['wetting']}}, -# 5: {0: {'wetting': exact_solution[5]['wetting']}} -# } - write_to_file = { 'meshes_and_markers': True, 'L_iterations': True } # initialise LDD simulation class -simulation = ldd.LDDsimulation(tol=1E-14, debug=False, LDDsolver_tol=1E-7) -simulation.set_parameters(output_dir="./output/with_cutoff_function", +simulation = ldd.LDDsimulation(tol=1E-14, debug=debugflag, LDDsolver_tol=solver_tol) +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, @@ -670,12 +557,12 @@ simulation.set_parameters(output_dir="./output/with_cutoff_function", dirichletBC_expression_strings=dirichletBC, exact_solution=exact_solution, densities=densities, - include_gravity=True, + include_gravity=include_gravity, write2file=write_to_file, ) simulation.initialise() # print(simulation.__dict__) -simulation.run() +simulation.run(analyse_condition=analyse_condition) # simulation.LDDsolver(time=0, debug=True, analyse_timestep=True) # df.info(parameters, True) diff --git a/TP-R-multi-patch-same-wetting-phase-as-RR-zero-nonwetting/TP-R-multi-patch-same-wetting-phase-as-RR-zero-nonwetting.py b/TP-R-multi-patch-same-wetting-phase-as-RR-zero-nonwetting/TP-R-multi-patch-same-wetting-phase-as-RR-zero-nonwetting.py index a808b2831b31460b2eb2b80b83af2046a9cf6b3c..80bde1a1ae5c45779b0dc4dc01a896d779307016 100755 --- a/TP-R-multi-patch-same-wetting-phase-as-RR-zero-nonwetting/TP-R-multi-patch-same-wetting-phase-as-RR-zero-nonwetting.py +++ b/TP-R-multi-patch-same-wetting-phase-as-RR-zero-nonwetting/TP-R-multi-patch-same-wetting-phase-as-RR-zero-nonwetting.py @@ -17,7 +17,7 @@ use_case = "TP-R-multi-patch-zero-nonwetting-phase" solver_tol = 1E-6 ############ GRID #######################ü -mesh_resolution = 50 +mesh_resolution = 40 timestep_size = 0.001 number_of_timesteps = 1500 # decide how many timesteps you want analysed. Analysed means, that we write out @@ -25,7 +25,7 @@ number_of_timesteps = 1500 number_of_timesteps_to_analyse = 10 starttime = 0 -Lw = 0.25 #/timestep_size +Lw = 0.5 #/timestep_size Lnw=Lw l_param_w = 40 diff --git a/TP-R-two-patch-test-case/TP-R-2-patch-test.py b/TP-R-two-patch-test-case/TP-R-2-patch-test.py index 168870d1569f1ce255675f364d4aa152b23a2592..62645f74fa57ed39e36b56632563868bb472c1ea 100755 --- a/TP-R-two-patch-test-case/TP-R-2-patch-test.py +++ b/TP-R-two-patch-test-case/TP-R-2-patch-test.py @@ -13,29 +13,29 @@ import helpers as hlp # init sympy session sym.init_printing() -use_case = "TP-R-two-patch" -solver_tol = 5E-7 +use_case = "TP-R-two-patch-realistic " +solver_tol = 1E-6 ############ GRID #######################ü -mesh_resolution = 40 -timestep_size = 0.000001 +mesh_resolution = 20 +timestep_size = 0.00001 number_of_timesteps = 15 # 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 = 5 starttime = 0 -Lw = 5 #/timestep_size +Lw = 1 #/timestep_size Lnw=Lw -l_param_w = 100 -l_param_nw = 100 +l_param_w = 40 +l_param_nw = 40 include_gravity = True debugflag = True analyse_condition = False -output_string = "./output/after_reimplementing_gravity_term_number_of_timesteps{}_".format(number_of_timesteps) +output_string = "./output/after_reimplementing_gravity_term-realistic_number_of_timesteps{}_".format(number_of_timesteps) ##### Domain and Interface #### # global simulation domain domain diff --git a/TP-one-patch/TP-one-patch-purely-postive-pc.py b/TP-one-patch/TP-one-patch-purely-postive-pc.py index 5638ef08f9064e249530647010ee086c653a363b..783c93b1d71c53749078fa44158b41baaad924dc 100755 --- a/TP-one-patch/TP-one-patch-purely-postive-pc.py +++ b/TP-one-patch/TP-one-patch-purely-postive-pc.py @@ -17,10 +17,10 @@ sym.init_printing() solver_tol = 1E-6 ############ GRID ####################### - -mesh_resolution = 40 -timestep_size = 0.001 -number_of_timesteps = 10 +use_case="TP-one-patch-purely-positive-pc" +mesh_resolution = 30 +timestep_size = 0.0005 +number_of_timesteps = 2500 # 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