fix smaller bugs, like numbering of timesteps

LDDsimulation/LDDsimulation.py

@@ -61,8 +61,9 @@ class LDDsimulation(object):
         # df.parameters['std_out_all_processes'] = False
 
 
+        # df.set_log_level(df.LogLevel.WARNING)
         df.set_log_level(df.LogLevel.INFO)
-        # df.set_log_level(df.LogLevel.INFO)
+        # df.set_log_level(df.LogLevel.DEBUG)
 
         # CRITICAL  = 50, // errors that may lead to data corruption and suchlike
         # ERROR     = 40, // things that go boom
@@ -70,7 +71,7 @@ class LDDsimulation(object):
         # INFO      = 20, // information of general interest
         # PROGRESS  = 16, // what's happening (broadly)
         # TRACE     = 13, // what's happening (in detail)
-        # DBG       = 10  // sundry
+        # DEBUG     = 10, // sundry
 
         hlp.print_once("Set default parameter...\n")
         self.restart_file = None
@@ -146,8 +147,8 @@ class LDDsimulation(object):
         # subdomain. This gets populated by self._init_solution_files()
         self.solution_file = None
         ########## SOLVER DEFINITION AND PARAMETERS ########
-        # print("\nLinear Algebra Backends:")
-        # df.list_linear_algebra_backends()
+        print("\nLinear Algebra Backends:")
+        df.list_linear_algebra_backends()
         # # print("\nLinear Solver Methods:")
         df.list_linear_solver_methods()
         print("\n")
@@ -156,12 +157,13 @@ class LDDsimulation(object):
         # # print("\nPeconditioners for Krylov Solvers:")
         df.list_krylov_solver_preconditioners()
 
-        # df.LinearSolver_default_parameters()
+
 
         df.info(df.parameters, True)
 
-        self.solver_type_is_Kryov = True
+        self.solver_type_is_Kryov = False
         ### Define the linear solver to be used.
+        if self.solver_type_is_Kryov:
             self.solver = 'bicgstab'#'superlu' #'gmres'#'bicgstab' # biconjugate gradient stabilized method
             self.preconditioner = 'jacobi' #'hypre_amg' 'default' #'ilu'#'hypre_amg' # incomplete LU factorization
             # dictionary of solver parametrs. This is passed to self._init_subdomains,
@@ -172,7 +174,19 @@ class LDDsimulation(object):
                 'absolute_tolerance': 1E-14,
                 'relative_tolerance': 1E-12,
                 'maximum_iterations': 1000,
-            'report': False
+                'report': True,
+                'monitor_convergence': True
+            }
+        else:
+            self.solver = 'superlu' #'gmres'#'bicgstab' # biconjugate gradient stabilized method
+            self.preconditioner = 'default' #'hypre_amg' 'default' #'ilu'#'hypre_amg' # incomplete LU factorization
+            # dictionary of solver parametrs. This is passed to self._init_subdomains,
+            # where for each subdomain a sovler object of type self.solver is created
+            # with these parameters.
+            self.solver_parameters = {
+                'report' : True,
+                'symmetric': False,
+                'verbose': True
             }
 
         self.debug_solver = debug
@@ -310,10 +324,10 @@ class LDDsimulation(object):
 
         df.info(colored("Start time stepping","green"))
         # write initial values to file.
-        self.timestep_num = int(0)
+        self.timestep_num = int(1)
         # 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:
+        while self.timestep_num <= self.number_of_timesteps:
             print(f"entering timestep ",
                   "{}".format(self.timestep_num),
                   "at time t = "+ "{number:.{digits}f}".format(number = self.t, digits = 4))
@@ -423,7 +437,7 @@ class LDDsimulation(object):
                                   )
                 # bypass subsequent error calculation for the first iteration, becauses
                 # it is zero anyways.
-                if iteration > 1:
+                #if iteration >= 1:
                 subsequent_iter_error = dict()
                 for phase in subdomain.has_phases:
                     error = df.Function(subdomain.function_space["pressure"][phase])
@@ -433,7 +447,7 @@ class LDDsimulation(object):
                         {phase: error_calculated}
                         )
                     if debug:
-                            print(f"time = {time}: subsequent error on subdomain {sd_index} for {phase} phase in iteration {iteration} = ", subsequent_iter_error[phase])
+                        print(f"time = {time} (at timestep {self.timestep_num}): subsequent error on subdomain {sd_index} for {phase} phase in iteration {iteration} = ", subsequent_iter_error[phase])
 
                 # calculate the maximum over phase of subsequent errors for the
                 # stopping criterion.
@@ -472,7 +486,7 @@ class LDDsimulation(object):
             # end loop over subdomains
             ##### stopping criterion for the solver.
             # only check if error criterion has been met after at least one iteration.
-            if iteration > 1:
+            #if iteration >= 1:
             if debug:
                 # print(f" total_subsequent_error in iter {iteration} = {total_subsequent_iter_error }\n")
                 print(f"the maximum distance between any to points in a cell mesh.hmax() on subdomain{sd_index} is h={subdomain.mesh.hmax()}")
@@ -966,16 +980,17 @@ class LDDsimulation(object):
             # df.LinearSolver(self.solver)
             if self.solver_type_is_Kryov:
                 self.subdomain[subdom_num].linear_solver = df.KrylovSolver(self.solver, self.preconditioner)
-                # we use the previous iteration as an initial guess for the linear solver.
-                solver_param = self.subdomain[subdom_num].linear_solver.parameters
-                df.info(solver_param, True)
-                solver_param['nonzero_initial_guess'] = self.solver_parameters['nonzero_initial_guess']
-                solver_param['absolute_tolerance'] = self.solver_parameters['absolute_tolerance']
-                solver_param['relative_tolerance'] = self.solver_parameters['relative_tolerance']
-                solver_param['maximum_iterations'] = self.solver_parameters['maximum_iterations']
-                solver_param['report'] = self.solver_parameters['report']
             else:
                 self.subdomain[subdom_num].linear_solver = df.LUSolver()
+                print("solver method:")
+                for parameter, value in self.subdomain[subdom_num].linear_solver.parameters.items():
+                    print(f"\'{parameter}\': {value}")
+
+            # assign solver parameters set at the beginning of the LDDsimulation class.
+            solver_param = self.subdomain[subdom_num].linear_solver.parameters
+            df.info(solver_param, True)
+            for parameter_name in self.solver_parameters.keys():
+                solver_param[parameter_name] = self.solver_parameters[parameter_name]
 
             # ## print out set solver parameters
             for parameter, value in self.subdomain[subdom_num].linear_solver.parameters.items():
@@ -1129,7 +1144,7 @@ class LDDsimulation(object):
             interpolation_degree = self.interpolation_degree
 
         if time is None:
-            time = self.starttime
+            time = self.t
         subdomain = self.subdomain[subdomain_index]
         # in case time == self.starttime, the expression has already been
         # assembled by self._init_DirichletBC_dictionary.
@@ -1215,7 +1230,7 @@ class LDDsimulation(object):
         if self.number_of_timesteps_to_analyse == 0:
             self.analyse_timesteps = None
         elif self.number_of_timesteps_to_analyse == 1:
-            self.analyse_timesteps = [0]
+            self.analyse_timesteps = [1]
         else:
             if not self.number_of_timesteps%(self.number_of_timesteps_to_analyse-1) == 0:
                 # if number_of_timesteps%(number_of_timesteps_to_analyse-1) is not 0, we
@@ -1225,7 +1240,7 @@ class LDDsimulation(object):
                 new_analysing_interval = int(round(self.number_of_timesteps/(self.number_of_timesteps_to_analyse-1)))
                 self.number_of_timesteps = new_analysing_interval*(self.number_of_timesteps_to_analyse)
                 print(f"I'm resetting number_of_timesteps = {self.number_of_timesteps}\n")
-            self.analyse_timesteps = np.linspace(0, self.number_of_timesteps, self.number_of_timesteps_to_analyse, dtype=int)
+            self.analyse_timesteps = np.linspace(1, self.number_of_timesteps, self.number_of_timesteps_to_analyse, dtype=int)
             print(f"the following timesteps will be analysed: {self.analyse_timesteps}")
 
         self.endtime = self.number_of_timesteps*self.timestep_size

TP-one-patch/TP-one-patch.py

@@ -17,9 +17,9 @@ sym.init_printing()
 solver_tol = 5E-6
 
 ############ GRID #######################ü
-mesh_resolution = 30
+mesh_resolution = 6
 timestep_size = 0.001
-number_of_timesteps = 600
+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
@@ -32,6 +32,7 @@ l_param_w = 80
 l_param_nw = 80
 
 include_gravity = True
+debugflag = False
 
 ##### Domain and Interface ####
 # global simulation domain domain
@@ -378,10 +379,10 @@ write_to_file = {
     'L_iterations': True
 }
 
-output_string = "./output/number_of_timesteps{}_".format(number_of_timesteps)
+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=False)
+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,#