diff --git a/.gitignore b/.gitignore index eb49951a7ee8b891af5d1b36525ca1d59b9c2ed5..d767444d7523f4672336d829b48089e41b21b0a7 100644 --- a/.gitignore +++ b/.gitignore @@ -60,6 +60,7 @@ *.ipynb core + # Ignoriere Bilder und Graphiken sowie Videos und Musik *.png *.jpg diff --git a/TP-TP-layered-soil-case-with-inner-patch-constant-solution/TP-TP-layered_soil_with_inner_patch_const_solution.py b/TP-TP-layered-soil-case-with-inner-patch-constant-solution/TP-TP-layered_soil_with_inner_patch_const_solution.py index d666eac6a424828bb4d120be87d496cf760684c5..a12790d75f3426f6b120d6dde4662605a0fa7eca 100755 --- a/TP-TP-layered-soil-case-with-inner-patch-constant-solution/TP-TP-layered_soil_with_inner_patch_const_solution.py +++ b/TP-TP-layered-soil-case-with-inner-patch-constant-solution/TP-TP-layered_soil_with_inner_patch_const_solution.py @@ -23,12 +23,12 @@ sym.init_printing() # ----------------------------------------------------------------------------# # ------------------- MESH ---------------------------------------------------# # ----------------------------------------------------------------------------# -mesh_resolution = 14 +mesh_resolution = 40 # ----------------------------------------:-----------------------------------# # ------------------- TIME ---------------------------------------------------# # ----------------------------------------------------------------------------# timestep_size = 0.0001 -number_of_timesteps = 10 +number_of_timesteps = 100 # 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 diff --git a/TP-TP-layered-soil-case/TP-TP-layered_soil-second-example.py b/TP-TP-layered-soil-case/TP-TP-layered_soil-second-example.py new file mode 100755 index 0000000000000000000000000000000000000000..23d75d3f8d88c427e7d174f73493e4f930d0c52d --- /dev/null +++ b/TP-TP-layered-soil-case/TP-TP-layered_soil-second-example.py @@ -0,0 +1,537 @@ +#!/usr/bin/python3 +"""This program sets up a domain together with a decomposition into subdomains +modelling layered soil. This is used for our LDD article with tp-tp and tp-r +coupling. + +Along with the subdomains and the mesh domain markers are set upself. +The resulting mesh is saved into files for later use. +""" + +#!/usr/bin/python3 +import dolfin as df +import mshr +import numpy as np +import sympy as sym +import typing as tp +import functools as ft +import domainPatch as dp +import LDDsimulation as ldd +import helpers as hlp + +# init sympy session +sym.init_printing() + +use_case="TP-TP-layered-soil-2nd-example" +solver_tol = 1E-5 + +############ GRID #######################ü +mesh_resolution = 5 +timestep_size = 0.001 +number_of_timesteps = 10 +# 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 = 0 +starttime = 0 + +Lw = 0.25 #/timestep_size +Lnw=Lw + +l_param_w = 40 +l_param_nw = 40 + +include_gravity = True +debugflag = False +analyse_condition = True + +output_string = "./output/test_postprocessing_number_of_timesteps{}_".format(number_of_timesteps) + +# global domain +subdomain0_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)] + +interface12_vertices = [df.Point(-1.0, 0.8), + df.Point(0.3, 0.8), + df.Point(0.5, 0.9), + df.Point(0.8, 0.7), + df.Point(1.0, 0.65)] +# subdomain1. +subdomain1_vertices = [interface12_vertices[0], + interface12_vertices[1], + interface12_vertices[2], + interface12_vertices[3], + interface12_vertices[4], # southern boundary, 12 interface + subdomain0_vertices[2], # eastern boundary, outer boundary + subdomain0_vertices[3]] # northern boundary, outer on_boundary + +# 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[4], # + subdomain0_vertices[2], # eastern boundary, outer boundary + subdomain0_vertices[3], + interface12_vertices[0]] +} + + +# interface23 +interface23_vertices = [df.Point(-1.0, 0.0), + df.Point(-0.35, 0.0), + # df.Point(6.5, 4.5), + df.Point(0.0, 0.0), + df.Point(0.5, 0.0), + # df.Point(11.5, 3.5), + # df.Point(13.0, 3) + df.Point(0.85, 0.0), + df.Point(1.0, 0.0) + ] + +#subdomain1 +subdomain2_vertices = [interface23_vertices[0], + interface23_vertices[1], + interface23_vertices[2], + interface23_vertices[3], + interface23_vertices[4], + interface23_vertices[5], # southern boundary, 23 interface + subdomain1_vertices[4], # eastern boundary, outer boundary + subdomain1_vertices[3], + subdomain1_vertices[2], + subdomain1_vertices[1], + subdomain1_vertices[0] ] # northern boundary, 12 interface + +subdomain2_outer_boundary_verts = { + 0: [interface23_vertices[5], + subdomain1_vertices[4]], + 1: [subdomain1_vertices[0], + interface23_vertices[0]] +} + + +# interface34 +interface34_vertices = [df.Point(-1.0, -0.6), + df.Point(-0.6, -0.45), + df.Point(0.3, -0.25), + df.Point(0.65, -0.6), + df.Point(1.0, -0.7)] + +# subdomain3 +subdomain3_vertices = [interface34_vertices[0], + interface34_vertices[1], + interface34_vertices[2], + interface34_vertices[3], + interface34_vertices[4], # southern boundary, 34 interface + subdomain2_vertices[5], # eastern boundary, outer boundary + subdomain2_vertices[4], + subdomain2_vertices[3], + subdomain2_vertices[2], + subdomain2_vertices[1], + subdomain2_vertices[0] ] # northern boundary, 23 interface + +subdomain3_outer_boundary_verts = { + 0: [interface34_vertices[4], + subdomain2_vertices[5]], + 1: [subdomain2_vertices[0], + interface34_vertices[0]] +} + +# subdomain4 +subdomain4_vertices = [subdomain0_vertices[0], + subdomain0_vertices[1], # southern boundary, outer boundary + subdomain3_vertices[4],# eastern boundary, outer boundary + subdomain3_vertices[3], + subdomain3_vertices[2], + subdomain3_vertices[1], + subdomain3_vertices[0] ] # northern boundary, 34 interface + +subdomain4_outer_boundary_verts = { + 0: [subdomain4_vertices[6], + subdomain4_vertices[0], + subdomain4_vertices[1], + subdomain4_vertices[2]] +} + + +subdomain_def_points = [subdomain0_vertices,# + subdomain1_vertices,# + subdomain2_vertices,# + subdomain3_vertices,# + subdomain4_vertices + ] + + +# interface_vertices introduces a global numbering of interfaces. +interface_def_points = [interface12_vertices, interface23_vertices, interface34_vertices] +adjacent_subdomains = [[1,2], [2,3], [3,4]] + +# 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, + 3: subdomain3_outer_boundary_verts, + 4: subdomain4_outer_boundary_verts +} + +isRichards = { + 1: False, + 2: False, + 3: False, + 4: False + } + +# isRichards = { +# 1: True, +# 2: True, +# 3: True, +# 4: True +# } + +# Dict of the form: { subdom_num : viscosity } +viscosity = { + 1: {'wetting' :1, + 'nonwetting': 1/50}, + 2: {'wetting' :1, + 'nonwetting': 1/50}, + 3: {'wetting' :1, + 'nonwetting': 1/50}, + 4: {'wetting' :1, + 'nonwetting': 1/50}, +} + +# Dict of the form: { subdom_num : density } +densities = { + 1: {'wetting': 997, #997 + 'nonwetting':1.225}, #1.225}}, + 2: {'wetting': 997, #997 + 'nonwetting':1.225}, #1.225}}, + 3: {'wetting': 997, #997 + 'nonwetting':1.225}, #1.225}}, + 4: {'wetting': 997, #997 + 'nonwetting':1.225}, #1.225}} +} + +gravity_acceleration = 9.81 +# porosities taken from +# https://www.geotechdata.info/parameter/soil-porosity.html +# Dict of the form: { subdom_num : porosity } +porosity = { + 1: 0.2, #0.2, # Clayey gravels, clayey sandy gravels + 2: 0.22, #0.22, # Silty gravels, silty sandy gravels + 3: 0.37, #0.37, # Clayey sands + 4: 0.2, #0.2 # Silty or sandy clay +} + +# subdom_num : subdomain L for L-scheme +L = { + 1: {'wetting' :Lw, + 'nonwetting': Lnw}, + 2: {'wetting' :Lw, + 'nonwetting': Lnw}, + 3: {'wetting' :Lw, + 'nonwetting': Lnw}, + 4: {'wetting' :Lw, + 'nonwetting': Lnw} +} + +# subdom_num : lambda parameter for the L-scheme +lambda_param = { + 1: {'wetting': l_param_w, + 'nonwetting': l_param_nw},# + 2: {'wetting': l_param_w, + 'nonwetting': l_param_nw},# + 3: {'wetting': l_param_w, + 'nonwetting': l_param_nw},# + 4: {'wetting': l_param_w, + 'nonwetting': l_param_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 + + +## 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_perm1w = ft.partial(rel_perm1w) +_rel_perm1nw = ft.partial(rel_perm1nw) +_rel_perm2w = ft.partial(rel_perm2w) +_rel_perm2nw = ft.partial(rel_perm2nw) + +subdomain1_rel_perm = { + 'wetting': _rel_perm1w,# + 'nonwetting': _rel_perm1nw +} + +subdomain2_rel_perm = { + 'wetting': _rel_perm2w,# + 'nonwetting': _rel_perm2nw +} + +# _rel_perm3 = ft.partial(rel_perm2) +# subdomain3_rel_perm = subdomain2_rel_perm.copy() +# +# _rel_perm4 = ft.partial(rel_perm1) +# subdomain4_rel_perm = subdomain1_rel_perm.copy() + +# dictionary of relative permeabilties on all domains. +relative_permeability = { + 1: subdomain1_rel_perm, + 2: subdomain1_rel_perm, + 3: subdomain2_rel_perm, + 4: 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: subdomain1_rel_perm_prime, + 3: subdomain2_rel_perm_prime, + 4: subdomain2_rel_perm_prime +} + + + +# 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 +# this function needs to be monotonically decreasing in the capillary pressure pc. +# since in the richards case pc=-pw, this becomes as a function of pw a mono +# tonically INCREASING function like in our Richards-Richards paper. However +# since we unify the treatment in the code for Richards and two-phase, we need +# the same requierment +# for both cases, two-phase and Richards. +def saturation(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 saturation_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 saturation_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(saturation_sym, n_index=3, alpha=0.001), + 2: ft.partial(saturation_sym, n_index=3, alpha=0.001), + 3: ft.partial(saturation_sym, n_index=6, alpha=0.001), + 4: ft.partial(saturation_sym, n_index=6, alpha=0.001) +} + +S_pc_sym_prime = { + 1: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001), + 2: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001), + 3: ft.partial(saturation_sym_prime, n_index=6, alpha=0.001), + 4: ft.partial(saturation_sym_prime, n_index=6, alpha=0.001) +} + +sat_pressure_relationship = { + 1: ft.partial(saturation, n_index=3, alpha=0.001), + 2: ft.partial(saturation, n_index=3, alpha=0.001), + 3: ft.partial(saturation, n_index=6, alpha=0.001), + 4: ft.partial(saturation, n_index=6, alpha=0.001) +} + + +############################################# +# 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_2patch = { + 1: {'wetting': -5 - (1+t*t)*(1 + x*x + y*y), + 'nonwetting': -1-t*(1.1+y + x**2)}, # - sym.sqrt(2+t**2)*(1-y)**2}, + 2: {'wetting': -5.0 - (1.0 + t*t)*(1.0 + x*x), + 'nonwetting': -1-t*(1.1 + x**2)},# - sym.sqrt(2+t**2)*(1-y)**2}, +} + +p_e_sym = { + 1: {'wetting': p_e_sym_2patch[1]['wetting'], + 'nonwetting': p_e_sym_2patch[1]['nonwetting']}, + 2: {'wetting': p_e_sym_2patch[1]['wetting'], + 'nonwetting': p_e_sym_2patch[1]['nonwetting']}, + 3: {'wetting': p_e_sym_2patch[2]['wetting'], + 'nonwetting': p_e_sym_2patch[2]['nonwetting']}, + 4: {'wetting': p_e_sym_2patch[2]['wetting'], + 'nonwetting': p_e_sym_2patch[2]['nonwetting']} +} + + +# p_e_sym = { +# 1: {'wetting': 1.0 - (1.0 + t*t)*(10.0 + x*x + (y-5.0)*(y-5.0)), +# 'nonwetting': - 2 - t*(1 + (y-5.0) + x**2)**2 -sym.sqrt(2+t**2)*(1 + (y-5.0)) }, +# 2: {'wetting': 1.0 - (1.0 + t*t)*(10.0 + x*x + (y-5.0)*(y-5.0)), +# 'nonwetting': - 2 - t*(1 + (y-5.0) + x**2)**2 -sym.sqrt(2+t**2)*(1 + (y-5.0))}, +# 3: {'wetting': 1.0 - (1.0 + t*t)*(10.0 + x*x + (y-5.0)*(y-5.0)) - (y-5.0)*(y-5.0)*3*sym.sin(-2*t+2*x)*sym.sin(1/2*y-1.2*t), +# 'nonwetting': - 2 - t*(1 + x**2)**2 -sym.sqrt(2+t**2)}, +# 4: {'wetting': 1.0 - (1.0 + t*t)*(10.0 + x*x + (y-5.0)*(y-5.0)) - (y-5.0)*(y-5.0)*3*sym.sin(-2*t+2*x)*sym.sin(1/2*y-1.2*t), +# 'nonwetting': - 2 - t*(1 + x**2)**2 -sym.sqrt(2+t**2)} +# } + +pc_e_sym = dict() +for subdomain, isR in isRichards.items(): + if isR: + pc_e_sym.update({subdomain: -p_e_sym[subdomain]['wetting']}) + else: + 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() +# 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]} + ) + +write_to_file = { + 'meshes_and_markers': True, + 'L_iterations': True +} + +# initialise LDD simulation class +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, + 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() +# 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/TP-TP-layered-soil-case/TP-TP-layered_soil.py b/TP-TP-layered-soil-case/TP-TP-layered_soil.py index 68d799f38cc340bda5969f610539b93b6d682455..a56a37a439dc7890dab0ed24fe9c62c59f986bed 100755 --- a/TP-TP-layered-soil-case/TP-TP-layered_soil.py +++ b/TP-TP-layered-soil-case/TP-TP-layered_soil.py @@ -16,27 +16,35 @@ import typing as tp import functools as ft import domainPatch as dp import LDDsimulation as ldd +import helpers as hlp # init sympy session sym.init_printing() -# ----------------------------------------------------------------------------# -# ------------------- MESH ---------------------------------------------------# -# ----------------------------------------------------------------------------# -mesh_resolution = 19 -# ----------------------------------------:-------------------------------------# -# ------------------- TIME ---------------------------------------------------# -# ----------------------------------------------------------------------------# -timestep_size = 0.001 -number_of_timesteps = 1000 -# 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 +use_case="TP-TP-layered-soil" +solver_tol = 5E-7 + +############ GRID #######################ü +mesh_resolution = 30 +timestep_size = 0.0005 +number_of_timesteps = 20 +# 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 = 5 starttime = 0 +Lw = 0.25 #/timestep_size +Lnw=Lw + l_param_w = 40 l_param_nw = 40 +include_gravity = True +debugflag = True +analyse_condition = False + +output_string = "./output/number_of_timesteps{}_".format(number_of_timesteps) + # global domain subdomain0_vertices = [df.Point(-1.0,-1.0), # df.Point(1.0,-1.0),# @@ -220,14 +228,14 @@ porosity = { # subdom_num : subdomain L for L-scheme L = { - 1: {'wetting' :0.25, - 'nonwetting': 0.25}, - 2: {'wetting' :0.25, - 'nonwetting': 0.25}, - 3: {'wetting' :0.25, - 'nonwetting': 0.25}, - 4: {'wetting' :0.25, - 'nonwetting': 0.25} + 1: {'wetting' :Lw, + 'nonwetting': Lnw}, + 2: {'wetting' :Lw, + 'nonwetting': Lnw}, + 3: {'wetting' :Lw, + 'nonwetting': Lnw}, + 4: {'wetting' :Lw, + 'nonwetting': Lnw} } # subdom_num : lambda parameter for the L-scheme @@ -302,7 +310,7 @@ def rel_perm1w_prime(s): def rel_perm1nw_prime(s): # relative permeabilty on subdomain1 - return 2*(1-s) + return -2*(1-s) # definition of the derivatives of the relative permeabilities # relative permeabilty functions on subdomain 1 @@ -312,7 +320,7 @@ def rel_perm2w_prime(s): def rel_perm2nw_prime(s): # relative permeabilty on subdomain1 - return 2*(1-s) + return -2*(1-s) _rel_perm1w_prime = ft.partial(rel_perm1w_prime) _rel_perm1nw_prime = ft.partial(rel_perm1nw_prime) @@ -398,10 +406,10 @@ x, y = sym.symbols('x[0], x[1]') # needed by UFL t = sym.symbols('t', positive=True) p_e_sym_2patch = { - 1: {'wetting': -1 - (1+t*t)*(1 + x*x + y*y), - 'nonwetting': -t*(1-y - x**2)**2 - sym.sqrt(2+t**2)*(1-y)}, - 2: {'wetting': -1.0 - (1.0 + t*t)*(1.0 + x*x), - 'nonwetting': -t*(1- x**2)**2 - sym.sqrt(2+t**2)*(1-y)}, + 1: {'wetting': -3 - (1+t*t)*(1 + x*x + y*y), + 'nonwetting': -1-t*(1-y - x**2)**2 - sym.sqrt(2+t**2)*(1-y)**2}, + 2: {'wetting': -3.0 - (1.0 + t*t)*(1.0 + x*x), + 'nonwetting': -1-t*(1- x**2)**2 - sym.sqrt(2+t**2)*(1-y)**2}, } p_e_sym = { @@ -427,83 +435,38 @@ p_e_sym = { # 'nonwetting': - 2 - t*(1 + x**2)**2 -sym.sqrt(2+t**2)} # } -pc_e_sym = { - 1: p_e_sym[1]['nonwetting'] - p_e_sym[1]['wetting'], - 2: p_e_sym[2]['nonwetting'] - p_e_sym[2]['wetting'], - 3: p_e_sym[3]['nonwetting'] - p_e_sym[3]['wetting'], - 4: p_e_sym[4]['nonwetting'] - p_e_sym[4]['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() @@ -539,8 +502,9 @@ write_to_file = { } # initialise LDD simulation class -simulation = ldd.LDDsimulation(tol=1E-14, debug=True, LDDsolver_tol=1E-7) -simulation.set_parameters(output_dir="./output/", +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, @@ -562,12 +526,12 @@ simulation.set_parameters(output_dir="./output/", 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)