SimpleNewtonNonlinearSolver.cpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2009
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(14);
00041     mTestDampingValues.push_back(-0.1);
00042     mTestDampingValues.push_back(0.05);
00043     for (unsigned i=1;i<=12;i++)
00044     {
00045         double val = double(i)/10;
00046         mTestDampingValues.push_back(val);
00047     }
00048 }
00049 
00050 SimpleNewtonNonlinearSolver::~SimpleNewtonNonlinearSolver()
00051 {}
00052 
00053 Vec SimpleNewtonNonlinearSolver::Solve(PetscErrorCode (*pComputeResidual)(SNES,Vec,Vec,void*),
00054                                        PetscErrorCode (*pComputeJacobian)(SNES,Vec,Mat*,Mat*,MatStructure*,void*),
00055                                        Vec initialGuess,
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);
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 #if (PETSC_VERSION_MINOR == 2) //Old API
00106             // VecWAXPY(&a, x, y, w) computes w = ax + y
00107             PetscScalar alpha = -mTestDampingValues[i];
00108             VecWAXPY(&alpha, negative_update, current_solution, test_vec);
00109 #else
00110             //[note: VecWAXPY(w,a,x,y) computes w = ax+y]
00111             VecWAXPY(test_vec,-mTestDampingValues[i],negative_update,current_solution);
00112 #endif
00113 
00114             // compute new residual
00115             linear_system.ZeroLinearSystem();
00116             (*pComputeResidual)(NULL, test_vec, linear_system.rGetRhsVector(), pContext);
00117             VecNorm(linear_system.rGetRhsVector(), NORM_2, &residual_norm);
00118             scaled_residual_norm = residual_norm/size;
00119 
00120             if (scaled_residual_norm < best_scaled_residual)
00121             {
00122                 best_scaled_residual = scaled_residual_norm;
00123                 best_damping_factor = mTestDampingValues[i];
00124             }
00125         }
00126         VecDestroy(test_vec);
00127 
00128         // check the smallest residual was actually smaller than the
00129         // previous. If not, quit.
00130         if (best_scaled_residual > old_scaled_residual_norm)
00131         {
00132             // Free memory
00133             VecDestroy(current_solution);
00134             VecDestroy(negative_update);
00135             // Raise error
00136             std::stringstream error_message;
00137             error_message << "Iteration " << counter << ", unable to find damping factor such that residual decreases in update direction";
00138             EXCEPTION(error_message.str());
00139         }
00140 
00141         if (mWriteStats)
00142         {
00143             std::cout << "    Best damping factor = " << best_damping_factor << "\n";
00144         }
00145 
00146 
00147         // update solution: current_guess = current_solution - best_damping_factor*negative_update
00148 #if (PETSC_VERSION_MINOR == 2) //Old API
00149         double minus_test_value = -best_damping_factor;
00150         VecAXPY(&minus_test_value, negative_update, current_solution);
00151 #else
00152         //[note: VecAXPY(y,a,x) computes y = ax+y]
00153         VecAXPY(current_solution, -best_damping_factor, negative_update);
00154 #endif
00155         scaled_residual_norm = best_scaled_residual;
00156         VecDestroy(negative_update);
00157 
00158         // compute best residual vector again and store in linear_system for next Solve()
00159         linear_system.ZeroLinearSystem();
00160         (*pComputeResidual)(NULL, current_solution, linear_system.rGetRhsVector(), pContext);
00161 
00162         if (mWriteStats)
00163         {
00164             std::cout << "    Iteration " << counter << ": ||residual||/N = " << scaled_residual_norm << "\n";
00165         }
00166     }
00167 
00168 
00169     if (mWriteStats)
00170     {
00171         std::cout << "  ..solved!\n\n";
00172     }
00173 
00174     return current_solution;
00175 }
00176 
00177 
00178 void SimpleNewtonNonlinearSolver::SetTolerance(double tolerance)
00179 {
00180     assert(tolerance > 0);
00181     mTolerance = tolerance;
00182 }
00183 

Generated on Wed Mar 18 12:51:52 2009 for Chaste by  doxygen 1.5.5