Fix errors related to the parameters access & setting. Update the tests.

src/bvpy/bvp.py

@@ -387,7 +387,15 @@ class BVP(object):
         parameters = self.solver.parameters._param
         krylov_parameter = self.solver.parameters._param['krylov_solver']
         for key in kwargs.keys():
-            if key in parameters:
+            if (key in krylov_parameter) and (key in parameters) and isinstance(self.solver, LinearSolver):
+                krylov_parameter[key] = kwargs[key]
+            elif key == 'krylov_solver':
+                for subkey in kwargs[key].keys():
+                    if subkey in krylov_parameter:
+                        krylov_parameter[subkey] = kwargs[key][subkey]
+                    else:
+                        logger.warn(str(subkey) + ' is not in the list of parameter for the krylov solver')
+            elif key in parameters:
                 parameters[key] = kwargs[key]
             elif key in krylov_parameter:
                 krylov_parameter[key] = kwargs[key]

src/bvpy/solvers/nonlinear.py

@@ -18,6 +18,7 @@
 #           https://www.gnu.org/licenses/lgpl-3.0.en.html
 #
 # -----------------------------------------------------------------------
+import logging
 
 import fenics as fe
 import time
@@ -180,16 +181,33 @@ class NonLinearSolver():
                 self._F, self._Jac, self._bnd, self._ps)
             self._solver = fe.PETScSNESSolver(MPI_COMM_WORLD)
 
-            paramters_dict = {}
-            paramters_dict.update(self.parameters['krylov_solver'])
-            paramters_dict['linear_solver'] = self.parameters['linear_solver']
+            for key, val in self.parameters._param.items():
+                if isinstance(val, dict):
+                    for subkey, subval in val.items():
+                        try:
+                            self._solver.parameters[key][subkey] = subval
+                        except:
+                            logging.warning(f'The parameter {key} {subkey} was not set!')
+                else:
+                    try:
+                        self._solver.parameters[key] = val
+                    except:
+                        logging.warning(f'The parameter {key} was not set!')
+
+            #paramters_dict = {}
+            #paramters_dict.update(self.parameters['krylov_solver'])
+            #paramters_dict['linear_solver'] = self.parameters['linear_solver']
             # TODO :: Check with Flo : is this normal that the NL solver has a 'linear_solver' param ?
-            paramters_dict['preconditioner'] = self.parameters['preconditioner']
-            paramters_dict['line_search'] = self.parameters['line_search']
-            paramters_dict['maximum_iterations'] = self.parameters['maximum_iterations']
-            monitor = paramters_dict.pop('monitor_convergence')
-            self._solver.parameters.update(paramters_dict)
-            self._solver.parameters['report'] = monitor
+            #paramters_dict['preconditioner'] = self.parameters['preconditioner']
+            #paramters_dict['line_search'] = self.parameters['line_search']
+            #paramters_dict['maximum_iterations'] = self.parameters['maximum_iterations']
+
+
+            #self._solver.parameters.update(paramters_dict)
+            #self._solver.parameters['report'] = monitor
+            #self._solver.parameters['krylov_solver'].update(self.parameters['krylov_solver']) # set the krylov parameters
+            #self._solver.parameters['absolute_tolerance']
+
             self._snes = self._solver.snes()
             self._snes.setConvergenceHistory()
 
@@ -203,17 +221,16 @@ class NonLinearSolver():
             of the performed computation.
         """
         if self._type == 'petsc':
-            if MPI_COMM_WORLD.Get_rank() == 0: # print only process of rank 0
-                logger.info(' Backend       : petsc')
-                logger.info(' Linear Solver : '
-                            + str(self._solver.parameters['linear_solver']))
-                logger.info(' Preconditioner: '
-                            + str(self._solver.parameters['preconditioner']))
-                logger.info(f" atol: {self._solver.parameters['absolute_tolerance']}")
-                logger.info(f" rtol: {self._solver.parameters['relative_tolerance']}")
-                logger.info(' Size          : '
-                            + str(self._u.function_space().dim()))
-                logger.info(' Solving NonLinearProblem ...')
+            logger.info(' Backend       : petsc')
+            logger.info(' Linear Solver : '
+                        + str(self._solver.parameters['linear_solver']))
+            logger.info(' Preconditioner: '
+                        + str(self._solver.parameters['preconditioner']))
+            logger.info(f" atol: {self._solver.parameters['absolute_tolerance']}")
+            logger.info(f" rtol: {self._solver.parameters['relative_tolerance']}")
+            logger.info(' Size          : '
+                        + str(self._u.function_space().dim()))
+            logger.info(' Solving NonLinearProblem ...')
 
             start = time.time()
             self._solver.solve(self._problem, self._u.vector())
@@ -221,10 +238,9 @@ class NonLinearSolver():
 
             residuals = self._snes.getConvergenceHistory()[0]
 
-            if MPI_COMM_WORLD.Get_rank() == 0:  # print only process of rank 0
-                logger.info(' ... Done [' + str(end - start) + ' s]')
-                conv = self._snes.getConvergedReason()
-                logger.info(' Iterations    : ' + str(len(residuals)))
-                logger.info(' Resiudal      : ' + str(residuals[-1]))
+            logger.info(' ... Done [' + str(end - start) + ' s]')
+            conv = self._snes.getConvergedReason()
+            logger.info(' Iterations    : ' + str(len(residuals)))
+            logger.info(' Resiudal      : ' + str(residuals[-1]))
 
             return residuals
\ No newline at end of file

src/bvpy/solvers/parameter.py

@@ -63,7 +63,8 @@ class SolverParameter():
                         'maximum_iterations': 50,
                         'maximum_residual_evaluations': 2000,
                         'method': 'default',
-                        'preconditioner': 'default',
+                        'report': False,
+                        'preconditioner': 'none',
                         'relative_tolerance': 1e-9,
                         'sign': 'default',
                         'solution_tolerance': 0.000000000000}

test/test_elasticity.py

@@ -42,7 +42,8 @@ class TestLinearElasticityProblem(unittest.TestCase):
 
         diri = dirichlet([0, 0], boundary='near(x, 0)')
         bvp = BVP(self.domain, vform, diri)
-        bvp.solve()
+        bvp.solve(krylov_solver={'absolute_tolerance':1e-10}) # assert that we can set absolute_tolerance for the krylov_solver
+        bvp.solve(absolute_tolerance=1e-10) # assert that the simplest way is also valid (autoset krylov parameter)
 
         comp_def = bvp.solution(self.L, self.W / 2)[1]
         err_rel = abs(self.theo_def_body + comp_def) / self.theo_def_body
@@ -56,7 +57,7 @@ class TestLinearElasticityProblem(unittest.TestCase):
         print(vform._parameters)
         diri = dirichlet([0, 0], boundary='near(x, 0)')
         bvp = BVP(self.domain, vform, diri)
-        bvp.solve(linear_solver='cg')
+        bvp.solve(linear_solver='superlu') # try another linear solver
 
         comp_def = bvp.solution(self.L, self.W / 2)[1]
         err_rel = abs(self.theo_def_body + comp_def) / self.theo_def_body
@@ -70,7 +71,7 @@ class TestLinearElasticityProblem(unittest.TestCase):
                                  plane_stress=True)
         diri = dirichlet([0, 0], boundary='near(x, 0)')
         bvp = BVP(self.domain, vform, diri)
-        bvp.solve(linear_solver='cg', absolute_tolerance=1e-7)
+        bvp.solve(linear_solver='cg', absolute_tolerance=1e-7) # try another linear solver and update absolute_tol
 
         comp_def = bvp.solution(self.L, self.W / 2)[1]
         err_rel = abs(self.theo_def_body + comp_def) / self.theo_def_body

test/test_transport.py

@@ -73,7 +73,7 @@ class TestTransport(unittest.TestCase):
 
         k1 = 9
         k2 = 11
-        dt = 0.1
+        dt = 0.05
 
         vform = CoupledTransportForm()
         vform.add_species('A', 1, lambda a, b: k1*(b-(a*b)/(1+b*b)))
@@ -104,4 +104,4 @@ class TestTransport(unittest.TestCase):
         ibvp.integrate(0, 1, dt,
                        linear_solver='mumps',
                        absolute_tolerance=1e-10,
-                       relative_tolerance=1e-10)
+                       relative_tolerance=1e-10, report=True)