From dc4c201ad3c2f9921fc67be1c3a421221db3fff5 Mon Sep 17 00:00:00 2001 From: David Seus <david.seus@ians.uni-stuttgart.de> Date: Wed, 7 Aug 2019 14:11:18 +0200 Subject: [PATCH] add analysis of condition number to code --- LDDsimulation/LDDsimulation.py | 91 +++++++++++++++---- TP-R-two-patch-test-case/TP-R-2-patch-test.py | 2 +- TP-one-patch/TP-one-patch.py | 41 +++++---- 3 files changed, 97 insertions(+), 37 deletions(-) diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py index b630fa5..e75dccb 100644 --- a/LDDsimulation/LDDsimulation.py +++ b/LDDsimulation/LDDsimulation.py @@ -5,6 +5,7 @@ LDD simulation class import os import dolfin as df import numpy as np +import numpy.linalg as la import typing as tp import mshr import domainPatch as dp @@ -293,7 +294,9 @@ class LDDsimulation(object): self._initialised = True - def run(self, analyse_solver_at_timesteps: tp.List[float] = None): + def run(self, + analyse_solver_at_timesteps: tp.List[float] = None, + analyse_condition: bool = False): """ run the simulation This method implements the time stepping, calls the acutal solvers and writes @@ -335,9 +338,13 @@ class LDDsimulation(object): # check if the solver should beanalised or not. if analyse_solver_at_timesteps is not None and np.isin(self.timestep_num, analyse_solver_at_timesteps, assume_unique=True): print("analyising timestep ...") - self.LDDsolver(debug = self.debug_solver, analyse_timestep = True) + self.LDDsolver(debug=self.debug_solver, + analyse_timestep=True, + analyse_condition=analyse_condition) else: - self.LDDsolver(debug = self.debug_solver, analyse_timestep = False) + self.LDDsolver(debug=self.debug_solver, + analyse_timestep=False, + analyse_condition=False) self.t += self.timestep_size self.timestep_num += 1 @@ -391,7 +398,8 @@ class LDDsimulation(object): def LDDsolver(self, time: float = None, # debug: bool = False, # - analyse_timestep: bool = False): + analyse_timestep: bool = False, + analyse_condition: bool = False): """ calulate the solution on all patches for the next timestep The heart of the simulation. The function LDDsolver implements the LDD @@ -433,8 +441,10 @@ class LDDsimulation(object): # set subdomain iteration to current iteration subdomain.iteration_number = iteration # solve the problem on subdomain - self.Lsolver_step(subdomain_index = sd_index,# - debug = debug, + condition_nums = self.Lsolver_step( + subdomain_index=sd_index,# + debug=debug, + analyse_condition=analyse_condition ) # bypass subsequent error calculation for the first iteration, becauses # it is zero anyways. @@ -464,8 +474,20 @@ class LDDsimulation(object): self.write_subsequent_errors_to_csv( filename = subsequent_error_filename, # subdomain_index = sd_index, - errors = subsequent_iter_error, - time = time) + errors = subsequent_iter_error + ) + if analyse_condition: + # save the calculated condition numbers of the assembled + # matrices a separate file for monitoring + condition_number_filename = self.output_dir\ + +self.output_filename_parameter_part[sd_index]\ + +"condition_numbers" +"_at_timestep"+\ + "{number}".format(number=self.timestep_num) +".csv" + self.write_condition_numbers_to_csv( + filename = condition_number_filename, # + subdomain_index = sd_index, + condition_number = condition_nums + ) # prepare next iteration # write the newly calculated solution to the inteface dictionaries @@ -496,7 +518,7 @@ class LDDsimulation(object): total_subsequent_error = np.amax(max_subsequent_error) if debug: - print(f"time = {time}: max subsequent error of both subdomain in iteration {iteration} = ", total_subsequent_error) + print(f"time = {time} (at timestep {self.timestep_num}):: max subsequent error of both subdomain in iteration {iteration} = ", total_subsequent_error) if total_subsequent_error < error_criterion: # if debug: print(f"all phases on all subdomains reached a subsequent error = {total_subsequent_error} < tol = {error_criterion} after {iteration} iterations.") @@ -602,7 +624,8 @@ class LDDsimulation(object): def Lsolver_step(self, subdomain_index: int,# - debug: bool = False) -> None: + analyse_condition: bool = False, + debug: bool = False) -> tp.Dict[str, float]: """ L-scheme solver iteration step for an object of class subdomain Lsolver_step implements L-scheme solver iteration step for an subdomain @@ -619,7 +642,9 @@ class LDDsimulation(object): """ subdomain = self.subdomain[subdomain_index] iteration = subdomain.iteration_number - + condition_number = None + if analyse_condition: + condition_number = dict() for phase in subdomain.has_phases: # extract L-scheme form and rhs (without gli term) from subdomain. governing_problem = subdomain.governing_problem(phase = phase) @@ -629,6 +654,11 @@ class LDDsimulation(object): form_assembled = df.assemble(form) rhs_assembled = df.assemble(rhs) + if analyse_condition: + condition_number.update( + {phase: la.cond(form_assembled.array())} + ) + # apply outer Dirichlet boundary conditions if present. if subdomain.outer_boundary is not None: # apply boundary dirichlet conditions for all outer boundaries. @@ -650,6 +680,8 @@ class LDDsimulation(object): if debug and subdomain.mesh.num_cells() < 36: print("\npressure after solver:\n", subdomain.pressure[phase].vector().get_local()) + return condition_number + def _init_solution_files(self): """ set up solution files for saving the solution of the LDD scheme for @@ -714,12 +746,8 @@ class LDDsimulation(object): filename: str, # subdomain_index: int,# errors: tp.Dict[str, float],# - time: float = None) -> None: + ) -> None: """ write subsequent errors to csv file for plotting with pgfplots""" - # in case no time is given, use the time of the class. - subdomain = self.subdomain[subdomain_index] - if time is None: - time = self.t subdomain = self.subdomain[subdomain_index] iteration = subdomain.iteration_number @@ -744,6 +772,37 @@ class LDDsimulation(object): else: self.errors.to_csv(filename, sep='\t', encoding='utf-8', index=False) + + def write_condition_numbers_to_csv(self,# + filename: str, # + subdomain_index: int,# + condition_number: tp.Dict[str, float],# + ) -> None: + """ write subsequent errors to csv file for plotting with pgfplots""" + subdomain = self.subdomain[subdomain_index] + iteration = subdomain.iteration_number + + if subdomain.isRichards: + condition = pd.DataFrame({ "iteration": iteration, + #"th. initial mass": df.assemble(rho*df.dx), \ + #"energy": self.total_energy(), \ + "wetting": condition_number["wetting"]}, index=[iteration]) + + else: + condition = pd.DataFrame({"iteration":iteration, #\ + #"th. initial mass": df.assemble(rho*df.dx), \ + #"energy": self.total_energy(), \ + "wetting": condition_number["wetting"], + "nonwetting": condition_number["nonwetting"]}, index=[iteration]) + # if self._rank == 0: + # check if file exists + if os.path.isfile(filename) == True: + with open(filename, 'a') as f: + condition.to_csv(f, header=False, sep='\t', encoding='utf-8', index=False) + else: + condition.to_csv(filename, sep='\t', encoding='utf-8', index=False) + + def write_errornorms_to_csv(self,# filename: str, # subdomain_index: int,# 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 ef696e3..9924115 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 @@ -14,7 +14,7 @@ import helpers as hlp sym.init_printing() solver_tol = 5e-7 -############ GRID #######################ü +######################## GRID ####################### mesh_resolution = 30 timestep_size = 0.0001 number_of_timesteps = 50 diff --git a/TP-one-patch/TP-one-patch.py b/TP-one-patch/TP-one-patch.py index c66fc33..4df57fc 100755 --- a/TP-one-patch/TP-one-patch.py +++ b/TP-one-patch/TP-one-patch.py @@ -14,25 +14,28 @@ import helpers as hlp sym.init_printing() -solver_tol = 5E-6 +solver_tol = 1E-7 ############ GRID #######################ü -mesh_resolution = 6 -timestep_size = 0.001 +mesh_resolution = 20 +timestep_size = 0.0001 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 = 10 starttime = 0 -Lw = 1/timestep_size +Lw = 1 #/timestep_size Lnw=Lw -l_param_w = 80 -l_param_nw = 80 +l_param_w = 40 +l_param_nw = 40 include_gravity = False -debugflag = False +debugflag = True +analyse_condition = False + +output_string = "./output/nondirichlet_number_of_timesteps{}_".format(number_of_timesteps) ##### Domain and Interface #### # global simulation domain domain @@ -283,16 +286,16 @@ cutoff = gaussian/(gaussian + zero_on_shrinking) # is_inside = ((1-sym.Abs(x) > epsilon) & (1-sym.Abs(y) > epsilon)) # sym.Piecewise((0, is_inside)) -# p_e_sym = { -# 0: {'wetting': (-3 - (1+t*t)*(1 + x*x + y*y))*cutoff, -# 'nonwetting': (-1 -t*(1+y + x**2)**2)*cutoff}, -# } - p_e_sym = { - 0: {'wetting': -(sym.cos(2*t-x - 2*y)*sym.sin(3*(1+y)/2*sym.pi)*sym.sin(5*(1+x)/2*sym.pi))**2, - 'nonwetting': -6*(sym.cos(t-x -y)*sym.sin(3*(1+y)/2*sym.pi)*sym.sin(5*(1+x)/2*sym.pi))**2}, + 0: {'wetting': (-3 - (1+t*t)*(1 + x*x + y*y)), #*cutoff, + 'nonwetting': (-1 -t*(1+y + x**2)**2)}, #*cutoff}, } +# p_e_sym = { +# 0: {'wetting': -(sym.cos(2*t-x - 2*y)*sym.sin(3*(1+y)/2*sym.pi)*sym.sin(5*(1+x)/2*sym.pi))**2, +# 'nonwetting': -6*(sym.cos(t-x -y)*sym.sin(3*(1+y)/2*sym.pi)*sym.sin(5*(1+x)/2*sym.pi))**2}, +# } + print(f"\n\n\nsymbolic type is {type(p_e_sym[0]['wetting'])}\n\n\n") # # pw_sym_x*pw_sym_y @@ -310,10 +313,10 @@ print(f"\n\n\nsymbolic type is {type(p_e_sym[0]['wetting'])}\n\n\n") pc_e_sym = dict() for subdomain, isR in isRichards.items(): if isR: - pc_e_sym.update({subdomain: -p_e_sym[subdomain]['wetting'].copy()}) + pc_e_sym.update({subdomain: -p_e_sym[subdomain]['wetting']}) else: - pc_e_sym.update({subdomain: p_e_sym[subdomain]['nonwetting'].copy() - - p_e_sym[subdomain]['wetting'].copy()}) + pc_e_sym.update({subdomain: p_e_sym[subdomain]['nonwetting'] + - p_e_sym[subdomain]['wetting']}) symbols = {"x": x, "y": y, @@ -379,8 +382,6 @@ write_to_file = { 'L_iterations': True } -output_string = "./output/testing_number_of_timesteps{}_".format(number_of_timesteps) - # initialise LDD simulation class simulation = ldd.LDDsimulation(tol = 1E-14, LDDsolver_tol=solver_tol, debug=debugflag) simulation.set_parameters(output_dir = output_string,# @@ -411,4 +412,4 @@ simulation.set_parameters(output_dir = output_string,# simulation.initialise() # simulation.write_exact_solution_to_xdmf() -simulation.run() +simulation.run(analyse_condition=analyse_condition) -- GitLab