diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py index aaa0b0d403ae9006b33cc90ed8d4b28018897a03..7b9b6698cea83c4b215120bf666db5b2655e04bc 100644 --- a/LDDsimulation/LDDsimulation.py +++ b/LDDsimulation/LDDsimulation.py @@ -106,7 +106,7 @@ class LDDsimulation(object): ## Private variables # maximal number of L-iterations that the LDD solver uses. - self._max_iter_num = 1000 + self._max_iter_num = 500 # TODO rewrite this with regard to the mesh sizes # self.calc_tol = self.tol # list of timesteps that get analyed. Gets initiated by self._init_analyse_timesteps @@ -198,6 +198,7 @@ class LDDsimulation(object): # done or not. def set_parameters(self,# + use_case: str, output_dir: str,# subdomain_def_points: tp.List[tp.List[df.Point]],# # this is actually an array of bools @@ -228,6 +229,7 @@ class LDDsimulation(object): )-> None: """ set parameters of an instance of class LDDsimulation""" + self.use_case = use_case self.output_dir = output_dir self.subdomain_def_points = subdomain_def_points self.isRichards = isRichards @@ -332,6 +334,7 @@ class LDDsimulation(object): # note that at this point of the code self.t = self.starttime as set in # the beginning. while self.timestep_num <= self.number_of_timesteps: + print("LDD simulation use case: {}".format(self.use_case)) print(f"entering timestep ", "{}".format(self.timestep_num), "at time t = "+ "{number:.{digits}f}".format(number = self.t, digits = 4)) @@ -347,7 +350,6 @@ class LDDsimulation(object): analyse_condition=False) self.t += self.timestep_size - self.timestep_num += 1 # write solutions to files. for subdom_ind, subdomain in self.subdomain.items(): subdomain.write_solution_to_xdmf(# @@ -380,8 +382,9 @@ class LDDsimulation(object): errors = relative_L2_errornorm, ) + self.timestep_num += 1 - + print("LDD simulation use case: {}".format(self.use_case)) df.info("Finished") df.info("Elapsed time:"+str(timer.elapsed()[0])) timer.stop() @@ -450,6 +453,8 @@ class LDDsimulation(object): # it is zero anyways. #if iteration >= 1: subsequent_iter_error = dict() + if debug: + print("LDD simulation use case: {}".format(self.use_case)) for phase in subdomain.has_phases: error = df.Function(subdomain.function_space["pressure"][phase]) error.assign(subdomain.pressure[phase] - subdomain.pressure_prev_iter[phase]) 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 ad5199e8eb04958b242dcef78bd247657c4ce0fe..f0710f0d532f03c1d7cdda0f1505c04ae6b66e79 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,28 +13,29 @@ import helpers as hlp # init sympy session sym.init_printing() -solver_tol = 1E-7 +use_case = "TP-R-two-patch" +solver_tol = 5E-7 ############ GRID #######################ü mesh_resolution = 40 -timestep_size = 0.0001 -number_of_timesteps = 10 +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 starttime = 0 -Lw = 4 #/timestep_size +Lw = 2.5 #/timestep_size Lnw=Lw -l_param_w = 60 -l_param_nw = 60 +l_param_w = 40 +l_param_nw = 40 include_gravity = True debugflag = True analyse_condition = False -output_string = "./output/nondirichlet_number_of_timesteps{}_".format(number_of_timesteps) +output_string = "./output/before_reimplementing_gravity_term_number_of_timesteps{}_".format(number_of_timesteps) ##### Domain and Interface #### # global simulation domain domain @@ -109,20 +110,20 @@ viscosity = {# 1 : {'wetting' :1}, #'nonwetting': 1}, # 2 : {'wetting' :1, - 'nonwetting': 1} + 'nonwetting': 1/50} } porosity = {# # subdom_num : porosity - 1 : 1,#0.22,# - 2 : 1#0.022 + 1 : 0.22,# + 2 : 0.0022 } # Dict of the form: { subdom_num : density } densities = { - 1: {'wetting': 1}, - 2: {'wetting': 1, - 'nonwetting': 1}, + 1: {'wetting': 997}, + 2: {'wetting': 997, + 'nonwetting': 1.225}, } gravity_acceleration = 9.81 @@ -322,9 +323,9 @@ x, y = sym.symbols('x[0], x[1]') # needed by UFL t = sym.symbols('t', positive=True) p_e_sym = { - 1: {'wetting': (-3.0 - (1.0 + t*t)*(1.0 + x*x + y*y))}, #*(1-x)**2*(1+x)**2*(1-y)**2}, - 2: {'wetting': (-3.0 - (1.0 + t*t)*(1.0 + x*x)), #*(1-x)**2*(1+x)**2*(1+y)**2, - 'nonwetting': (-t**2*(1+y + x**2)**2 - sym.sqrt(2+t**4))*y**2}, #*(1-x)**2*(1+x)**2*(1+y)**2}, + 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) - sym.sqrt(2+t**4))*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) @@ -403,8 +404,9 @@ write_to_file = { # initialise LDD simulation class -simulation = ldd.LDDsimulation(tol = 1E-14, LDDsolver_tol=solver_tol, debug=True) -simulation.set_parameters(output_dir = output_string,# +simulation = ldd.LDDsimulation(tol = 1E-14, LDDsolver_tol=solver_tol, debug=debugflag) +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,# diff --git a/TP-TP-2-patch-pure-dd/TP-TP-2-patch-pure-dd.py b/TP-TP-2-patch-pure-dd/TP-TP-2-patch-pure-dd.py index c9b8d4a58c6549db59de9fa92d9b6034a4e576c1..1a9a2c1ebb78ae2bd02e22d1ac3e8e7b1446aa0e 100755 --- a/TP-TP-2-patch-pure-dd/TP-TP-2-patch-pure-dd.py +++ b/TP-TP-2-patch-pure-dd/TP-TP-2-patch-pure-dd.py @@ -13,26 +13,27 @@ import helpers as hlp # init sympy session sym.init_printing() -solver_tol = 5E-7 +use_case = "TP-TP-2-patch-pure-dd" +solver_tol = 1E-6 ############ GRID #######################ü -mesh_resolution = 50 +mesh_resolution = 30 timestep_size = 0.0001 -number_of_timesteps = 1000 +number_of_timesteps = 11000 # 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 = 1 #/timestep_size +Lw = 10 #/timestep_size Lnw=Lw -l_param_w = 40 -l_param_nw = 40 +l_param_w = 20 +l_param_nw = 20 include_gravity = True debugflag = False -analyse_condition = True +analyse_condition = False output_string = "./output/nondirichlet_number_of_timesteps{}_".format(number_of_timesteps) @@ -347,85 +348,93 @@ symbols = { "x": x, "y": y, "t": t} -epsilon_x_inner = 0.7 -epsilon_x_outer = 0.99 -epsilon_y_inner = epsilon_x_inner -epsilon_y_outer = epsilon_x_outer - -def mollifier(x, epsilon): - """ one d mollifier """ - out_expr = sym.exp(-1/(1-(x/epsilon)**2) + 1) - return out_expr - -mollifier_handle = ft.partial(mollifier, epsilon=epsilon_x_inner) - -pw_sym_x = sym.Piecewise( - (mollifier_handle(x), x**2 < epsilon_x_outer**2), - (0, True) -) -pw_sym_y = sym.Piecewise( - (mollifier_handle(y), y**2 < epsilon_y_outer**2), - (0, True) -) - -def mollifier2d(x, y, epsilon): - """ one d mollifier """ - out_expr = sym.exp(-1/(1-(x**2 + y**2)/epsilon**2) + 1) - return out_expr - -mollifier2d_handle = ft.partial(mollifier2d, epsilon=epsilon_x_outer) - -pw_sym2d_x = sym.Piecewise( - (mollifier2d_handle(x, y), x**2 + y**2 < epsilon_x_outer**2), - (0, True) -) - -zero_on_epsilon_shrinking_of_subdomain = sym.Piecewise( - (mollifier_handle(sym.sqrt(x**2 + y**2)+2*epsilon_x_inner), ((-2*epsilon_x_inner<sym.sqrt(x**2 + y**2)) & (sym.sqrt(x**2 + y**2)<-epsilon_x_inner))), - (0, ((-epsilon_x_inner<=sym.sqrt(x**2 + y**2)) & (sym.sqrt(x**2 + y**2)<=epsilon_x_inner))), - (mollifier_handle(sym.sqrt(x**2 + y**2)-2*epsilon_x_inner), ((epsilon_x_inner<sym.sqrt(x**2 + y**2)) & (sym.sqrt(x**2 + y**2)<2*epsilon_x_inner))), - (1, True), -) - -zero_on_epsilon_shrinking_of_subdomain_x = sym.Piecewise( - (mollifier_handle(x+2*epsilon_x_inner), ((-2*epsilon_x_inner<x) & (x<-epsilon_x_inner))), - (0, ((-epsilon_x_inner<=x) & (x<=epsilon_x_inner))), - (mollifier_handle(x-2*epsilon_x_inner), ((epsilon_x_inner<x) & (x<2*epsilon_x_inner))), - (1, True), -) - -zero_on_epsilon_shrinking_of_subdomain_y = sym.Piecewise( - (1, y<=-2*epsilon_x_inner), - (mollifier_handle(y+2*epsilon_x_inner), ((-2*epsilon_x_inner<y) & (y<-epsilon_x_inner))), - (0, ((-epsilon_x_inner<=y) & (y<=epsilon_x_inner))), - (mollifier_handle(y-2*epsilon_x_inner), ((epsilon_x_inner<y) & (y<2*epsilon_x_inner))), - (1, True), -) - -zero_on_shrinking = zero_on_epsilon_shrinking_of_subdomain #zero_on_epsilon_shrinking_of_subdomain_x + zero_on_epsilon_shrinking_of_subdomain_y -gaussian = pw_sym2d_x# pw_sym_y*pw_sym_x -cutoff = gaussian/(gaussian + zero_on_shrinking) - - -sat_sym = { - 1: 0.5 + 0.25*sym.sin(x-t)*sym.cos(y-t), - 2: 0.5 + 0.25*sym.sin(x-t)*sym.cos(y-t) - } - -Spc = { - 1: sym.Piecewise((pc_saturation_sym[1](sat_sym[1]), sat_sym[1] > 0), (pc_saturation_sym[1](sat_sym[1]), 1>=sat_sym[1]), (0, True)), - 2: sym.Piecewise((pc_saturation_sym[2](sat_sym[2]), sat_sym[2] > 0), (pc_saturation_sym[2](sat_sym[2]), 2>=sat_sym[2]), (0, True)) - } +# epsilon_x_inner = 0.7 +# epsilon_x_outer = 0.99 +# epsilon_y_inner = epsilon_x_inner +# epsilon_y_outer = epsilon_x_outer +# +# def mollifier(x, epsilon): +# """ one d mollifier """ +# out_expr = sym.exp(-1/(1-(x/epsilon)**2) + 1) +# return out_expr +# +# mollifier_handle = ft.partial(mollifier, epsilon=epsilon_x_inner) +# +# pw_sym_x = sym.Piecewise( +# (mollifier_handle(x), x**2 < epsilon_x_outer**2), +# (0, True) +# ) +# pw_sym_y = sym.Piecewise( +# (mollifier_handle(y), y**2 < epsilon_y_outer**2), +# (0, True) +# ) +# +# def mollifier2d(x, y, epsilon): +# """ one d mollifier """ +# out_expr = sym.exp(-1/(1-(x**2 + y**2)/epsilon**2) + 1) +# return out_expr +# +# mollifier2d_handle = ft.partial(mollifier2d, epsilon=epsilon_x_outer) +# +# pw_sym2d_x = sym.Piecewise( +# (mollifier2d_handle(x, y), x**2 + y**2 < epsilon_x_outer**2), +# (0, True) +# ) +# +# zero_on_epsilon_shrinking_of_subdomain = sym.Piecewise( +# (mollifier_handle(sym.sqrt(x**2 + y**2)+2*epsilon_x_inner), ((-2*epsilon_x_inner<sym.sqrt(x**2 + y**2)) & (sym.sqrt(x**2 + y**2)<-epsilon_x_inner))), +# (0, ((-epsilon_x_inner<=sym.sqrt(x**2 + y**2)) & (sym.sqrt(x**2 + y**2)<=epsilon_x_inner))), +# (mollifier_handle(sym.sqrt(x**2 + y**2)-2*epsilon_x_inner), ((epsilon_x_inner<sym.sqrt(x**2 + y**2)) & (sym.sqrt(x**2 + y**2)<2*epsilon_x_inner))), +# (1, True), +# ) +# +# zero_on_epsilon_shrinking_of_subdomain_x = sym.Piecewise( +# (mollifier_handle(x+2*epsilon_x_inner), ((-2*epsilon_x_inner<x) & (x<-epsilon_x_inner))), +# (0, ((-epsilon_x_inner<=x) & (x<=epsilon_x_inner))), +# (mollifier_handle(x-2*epsilon_x_inner), ((epsilon_x_inner<x) & (x<2*epsilon_x_inner))), +# (1, True), +# ) +# +# zero_on_epsilon_shrinking_of_subdomain_y = sym.Piecewise( +# (1, y<=-2*epsilon_x_inner), +# (mollifier_handle(y+2*epsilon_x_inner), ((-2*epsilon_x_inner<y) & (y<-epsilon_x_inner))), +# (0, ((-epsilon_x_inner<=y) & (y<=epsilon_x_inner))), +# (mollifier_handle(y-2*epsilon_x_inner), ((epsilon_x_inner<y) & (y<2*epsilon_x_inner))), +# (1, True), +# ) +# +# zero_on_shrinking = zero_on_epsilon_shrinking_of_subdomain #zero_on_epsilon_shrinking_of_subdomain_x + zero_on_epsilon_shrinking_of_subdomain_y +# gaussian = pw_sym2d_x# pw_sym_y*pw_sym_x +# cutoff = gaussian/(gaussian + zero_on_shrinking) +# +# +# sat_sym = { +# 1: 0.5 + 0.25*sym.sin(x-t)*sym.cos(y-t), +# 2: 0.5 + 0.25*sym.sin(x-t)*sym.cos(y-t) +# } +# +# Spc = { +# 1: sym.Piecewise((pc_saturation_sym[1](sat_sym[1]), sat_sym[1] > 0), (pc_saturation_sym[1](sat_sym[1]), 1>=sat_sym[1]), (0, True)), +# 2: sym.Piecewise((pc_saturation_sym[2](sat_sym[2]), sat_sym[2] > 0), (pc_saturation_sym[2](sat_sym[2]), 2>=sat_sym[2]), (0, True)) +# } +# +# p1w = (-1 - (1+t*t)*(1 + x*x + y*y))#*cutoff +# p2w = p1w +# p_e_sym = { +# 1: {'wetting': p1w, +# 'nonwetting': (p1w + Spc[1])}, #*cutoff}, +# 2: {'wetting': p2w, +# 'nonwetting': (p2w + Spc[2])}, #*cutoff}, +# } -p1w = (-1 - (1+t*t)*(1 + x*x + y*y))#*cutoff -p2w = p1w p_e_sym = { - 1: {'wetting': p1w, - 'nonwetting': (p1w + Spc[1])}, #*cutoff}, - 2: {'wetting': p2w, - 'nonwetting': (p2w + Spc[2])}, #*cutoff}, + 1: {'wetting': (-5 - (1+t*t)*(1 + x*x + y*y)), #*cutoff, + 'nonwetting': (-1 -t*(1.1+y + x**2))}, #*cutoff}, + 2: {'wetting': (-5 - (1+t*t)*(1 + x*x + y*y)), #*cutoff, + 'nonwetting': (-1 -t*(1.1+y + x**2))}, #*cutoff}, } + pc_e_sym = dict() for subdomain, isR in isRichards.items(): if isR: @@ -434,10 +443,7 @@ for subdomain, isR in isRichards.items(): pc_e_sym.update({subdomain: p_e_sym[subdomain]['nonwetting'] - p_e_sym[subdomain]['wetting']}) -# pc_e_sym = { -# 1: -1*p_e_sym[1]['wetting'], -# 2: -1*p_e_sym[2]['wetting'], -# } + exact_solution_example = hlp.generate_exact_solution_expressions( symbols=symbols, @@ -501,31 +507,32 @@ write_to_file = { # initialise LDD simulation class simulation = ldd.LDDsimulation(tol=1E-14, LDDsolver_tol=solver_tol, debug=debugflag) -simulation.set_parameters(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.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() # simulation.write_exact_solution_to_xdmf() diff --git a/TP-TP-layered-soil-case-with-inner-patch/TP-TP-layered_soil_with_inner_patch.py b/TP-TP-layered-soil-case-with-inner-patch/TP-TP-layered_soil_with_inner_patch.py index b1a38f81d538ed998f946023cb60c7a711fff34d..126c91e6fa3985b138332740dec5a1c6cda3178a 100755 --- a/TP-TP-layered-soil-case-with-inner-patch/TP-TP-layered_soil_with_inner_patch.py +++ b/TP-TP-layered-soil-case-with-inner-patch/TP-TP-layered_soil_with_inner_patch.py @@ -21,12 +21,13 @@ import helpers as hlp # init sympy session sym.init_printing() +use_case = "TP-TP-layered-soil-with-inner-patch" solver_tol = 5E-7 ############ GRID #######################ü mesh_resolution = 40 -timestep_size = 0.0005 -number_of_timesteps = 1000 +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 = 5 @@ -38,7 +39,7 @@ Lnw=Lw l_param_w = 40 l_param_nw = 40 -include_gravity = True +include_gravity = False debugflag = False analyse_condition = True @@ -701,7 +702,8 @@ write_to_file = { # initialise LDD simulation class simulation = ldd.LDDsimulation(tol=1E-14, debug=debugflag, LDDsolver_tol=solver_tol) -simulation.set_parameters(output_dir=output_string, +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, 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 1d77fdd71a1b7e58556975b54aa74fea02e611b1..a56a37a439dc7890dab0ed24fe9c62c59f986bed 100755 --- a/TP-TP-layered-soil-case/TP-TP-layered_soil.py +++ b/TP-TP-layered-soil-case/TP-TP-layered_soil.py @@ -21,12 +21,13 @@ import helpers as hlp # init sympy session sym.init_printing() +use_case="TP-TP-layered-soil" solver_tol = 5E-7 ############ GRID #######################ü -mesh_resolution = 5 -timestep_size = 0.0001 -number_of_timesteps = 2000 +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 @@ -35,14 +36,14 @@ starttime = 0 Lw = 0.25 #/timestep_size Lnw=Lw -l_param_w = 60 -l_param_nw = 60 +l_param_w = 40 +l_param_nw = 40 include_gravity = True -debugflag = False -analyse_condition = True +debugflag = True +analyse_condition = False -output_string = "./output/brief_docker_test_number_of_timesteps{}_".format(number_of_timesteps) +output_string = "./output/number_of_timesteps{}_".format(number_of_timesteps) # global domain subdomain0_vertices = [df.Point(-1.0,-1.0), # @@ -502,7 +503,8 @@ write_to_file = { # initialise LDD simulation class simulation = ldd.LDDsimulation(tol=1E-14, debug=debugflag, LDDsolver_tol=solver_tol) -simulation.set_parameters(output_dir=output_string, +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, diff --git a/TP-one-patch/TP-one-patch.py b/TP-one-patch/TP-one-patch.py index dbf4cf3e8dc83fb7a33ae7a0f70f0b0dfb074538..ab050acbb1a298500b74c419f6b1710983e4e90b 100755 --- a/TP-one-patch/TP-one-patch.py +++ b/TP-one-patch/TP-one-patch.py @@ -19,7 +19,7 @@ solver_tol = 5E-6 ############ GRID #######################ü mesh_resolution = 30 timestep_size = 0.0005 -number_of_timesteps = 2500 +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 = 5