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)