SimpleNewtonNonlinearSolver.cpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2011
00004 
00005 University of Oxford means the Chancellor, Masters and Scholars of the
00006 University of Oxford, having an administrative office at Wellington
00007 Square, Oxford OX1 2JD, UK.
00008 
00009 This file is part of Chaste.
00010 
00011 Chaste is free software: you can redistribute it and/or modify it
00012 under the terms of the GNU Lesser General Public License as published
00013 by the Free Software Foundation, either version 2.1 of the License, or
00014 (at your option) any later version.
00015 
00016 Chaste is distributed in the hope that it will be useful, but WITHOUT
00017 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
00019 License for more details. The offer of Chaste under the terms of the
00020 License is subject to the License being interpreted in accordance with
00021 English Law and subject to any action against the University of Oxford
00022 being under the jurisdiction of the English Courts.
00023 
00024 You should have received a copy of the GNU Lesser General Public License
00025 along with Chaste. If not, see <http://www.gnu.org/licenses/>.
00026 
00027 */
00028 
00029 
00030 #include "SimpleNewtonNonlinearSolver.hpp"
00031 #include <iostream>
00032 #include <cassert>
00033 #include "Exception.hpp"
00034 
00035 SimpleNewtonNonlinearSolver::SimpleNewtonNonlinearSolver(double linearSolverRelativeTolerance)
00036         : mLinearSolverRelativeTolerance(linearSolverRelativeTolerance),
00037         mTolerance(1e-5),
00038         mWriteStats(false)
00039 {
00040     mTestDampingValues.push_back(-0.1);
00041     mTestDampingValues.push_back(0.05);
00042     for (unsigned i=1; i<=12; i++)
00043     {
00044         double val = double(i)/10;
00045         mTestDampingValues.push_back(val);
00046     }
00047 }
00048 
00049 SimpleNewtonNonlinearSolver::~SimpleNewtonNonlinearSolver()
00050 {}
00051 
00052 Vec SimpleNewtonNonlinearSolver::Solve(PetscErrorCode (*pComputeResidual)(SNES,Vec,Vec,void*),
00053                                        PetscErrorCode (*pComputeJacobian)(SNES,Vec,Mat*,Mat*,MatStructure*,void*),
00054                                        Vec initialGuess,
00055                                        unsigned fill,
00056                                        void* pContext)
00057 {
00058     PetscInt size;
00059     VecGetSize(initialGuess, &size);
00060 
00061     Vec current_solution;
00062     VecDuplicate(initialGuess, &current_solution);
00063     VecCopy(initialGuess, current_solution);
00064 
00065     LinearSystem linear_system(current_solution, fill);
00066 
00067     (*pComputeResidual)(NULL, current_solution, linear_system.rGetRhsVector(), pContext);
00068 
00069 
00070     double residual_norm;
00071     VecNorm(linear_system.rGetRhsVector(), NORM_2, &residual_norm);
00072     double scaled_residual_norm = residual_norm/size;
00073 
00074     if (mWriteStats)
00075     {
00076         std::cout << "Newton's method:\n  Initial ||residual||/N = " << scaled_residual_norm
00077                   << "\n  Attempting to solve to tolerance " << mTolerance << "..\n";
00078     }
00079 
00080     double old_scaled_residual_norm;
00081     unsigned counter = 0;
00082     while (scaled_residual_norm > mTolerance)
00083     {
00084         counter++;
00085 
00086         // store the old norm to check with the new later
00087         old_scaled_residual_norm = scaled_residual_norm;
00088 
00089         // compute Jacobian and solve J dx = f for the (negative) update dx, (J the jacobian, f the residual)
00090         (*pComputeJacobian)(NULL, current_solution, &(linear_system.rGetLhsMatrix()), NULL, NULL, pContext);
00091 
00092         Vec negative_update = linear_system.Solve();
00093 
00094 
00095         Vec test_vec;
00096         VecDuplicate(initialGuess, &test_vec);
00097 
00098         double best_damping_factor = 1.0;
00099         double best_scaled_residual = 1e20; //large
00100 
00101         // loop over all the possible damping value and determine which gives
00102         // smallest residual
00103         for (unsigned i=0; i<mTestDampingValues.size(); i++)
00104         {
00105             //note: WAXPY calls VecWAXPY(w,a,x,y) which computes w = ax+y 
00106             PetscVecTools::WAXPY(test_vec,-mTestDampingValues[i],negative_update,current_solution);
00107 
00108             // compute new residual
00109             linear_system.ZeroLinearSystem();
00110             (*pComputeResidual)(NULL, test_vec, linear_system.rGetRhsVector(), pContext);
00111             VecNorm(linear_system.rGetRhsVector(), NORM_2, &residual_norm);
00112             scaled_residual_norm = residual_norm/size;
00113 
00114             if (scaled_residual_norm < best_scaled_residual)
00115             {
00116                 best_scaled_residual = scaled_residual_norm;
00117                 best_damping_factor = mTestDampingValues[i];
00118             }
00119         }
00120         VecDestroy(test_vec);
00121 
00122         // check the smallest residual was actually smaller than the
00123         // previous. If not, quit.
00124         if (best_scaled_residual > old_scaled_residual_norm)
00125         {
00126             // Free memory
00127             VecDestroy(current_solution);
00128             VecDestroy(negative_update);
00129             // Raise error
00130             std::stringstream error_message;
00131             error_message << "Iteration " << counter << ", unable to find damping factor such that residual decreases in update direction";
00132             EXCEPTION(error_message.str());
00133         }
00134 
00135         if (mWriteStats)
00136         {
00137             std::cout << "    Best damping factor = " << best_damping_factor << "\n";
00138         }
00139 
00140 
00141         // update solution: current_guess = current_solution - best_damping_factor*negative_update
00142         PetscVecTools::AddScaledVector(current_solution, negative_update, -best_damping_factor);
00143         scaled_residual_norm = best_scaled_residual;
00144         VecDestroy(negative_update);
00145 
00146         // compute best residual vector again and store in linear_system for next Solve()
00147         linear_system.ZeroLinearSystem();
00148         (*pComputeResidual)(NULL, current_solution, linear_system.rGetRhsVector(), pContext);
00149 
00150         if (mWriteStats)
00151         {
00152             std::cout << "    Iteration " << counter << ": ||residual||/N = " << scaled_residual_norm << "\n";
00153         }
00154     }
00155 
00156     if (mWriteStats)
00157     {
00158         std::cout << "  ..solved!\n\n";
00159     }
00160 
00161     return current_solution;
00162 }
00163 
00164 
00165 void SimpleNewtonNonlinearSolver::SetTolerance(double tolerance)
00166 {
00167     assert(tolerance > 0);
00168     mTolerance = tolerance;
00169 }
00170 

Generated on Tue May 31 14:31:42 2011 for Chaste by  doxygen 1.5.5