diff --git a/src/bvpy/bvp.py b/src/bvpy/bvp.py index b587a8299f00639e8a59991230cc1886ddab0e4a..0a9210e44c8c0d7db921df2a1aaa389ad5c6cc3e 100644 --- a/src/bvpy/bvp.py +++ b/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] diff --git a/src/bvpy/solvers/nonlinear.py b/src/bvpy/solvers/nonlinear.py index a4028bd4c4dd22a98627a3eb4051338c5d94294e..ed01e3f911b1aba629c90d94bb5ff59a02a6ee6b 100644 --- a/src/bvpy/solvers/nonlinear.py +++ b/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 diff --git a/src/bvpy/solvers/parameter.py b/src/bvpy/solvers/parameter.py index 3e50ca9043cd0572ef9030d97d2a2c1f3996e7c9..8f6af29060131250c615d9ae59236cc64c168609 100644 --- a/src/bvpy/solvers/parameter.py +++ b/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} diff --git a/test/test_elasticity.py b/test/test_elasticity.py index 8f2e8ca019b48090803e9e056f41db1a82de485c..20baecdf92903a14033ea25e87d05c1dff5828be 100644 --- a/test/test_elasticity.py +++ b/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 diff --git a/test/test_transport.py b/test/test_transport.py index 5f730ca8cc1da3e5d097178ebba8f82bdb0682be..c932ed64c22457cdd171a6f0bb45224d150d853f 100644 --- a/test/test_transport.py +++ b/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)