AbstractNonlinearAssemblerSolverHybrid.hpp

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 #ifndef ABSTRACTNONLINEARASSEMBLERSOLVERHYBRID_HPP_
00030 #define ABSTRACTNONLINEARASSEMBLERSOLVERHYBRID_HPP_
00031 
00032 
00033 
00034 #include "BoundaryConditionsContainer.hpp"
00035 #include "AbstractNonlinearEllipticPde.hpp"
00036 #include "AbstractNonlinearSolver.hpp"
00037 #include "PetscTools.hpp"
00038 #include "PetscMatTools.hpp"
00039 #include "PetscVecTools.hpp"
00040 #include "AbstractFeObjectAssembler.hpp"
00041 #include "LinearSystem.hpp"
00042 #include "SimplePetscNonlinearSolver.hpp"
00043 
00044 
00055 
00056 
00071 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00072 PetscErrorCode AbstractNonlinearAssemblerSolverHybrid_ComputeResidual(SNES snes,
00073                                                                       Vec currentGuess,
00074                                                                       Vec residualVector,
00075                                                                       void* pContext);
00076 
00077 
00078 
00095 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00096 PetscErrorCode AbstractNonlinearAssemblerSolverHybrid_ComputeJacobian(SNES snes,
00097                                                                       Vec currentGuess,
00098                                                                       Mat* pGlobalJacobian,
00099                                                                       Mat* pPreconditioner,
00100                                                                       MatStructure* pMatStructure,
00101                                                                       void* pContext);
00102 
00103 
00104 
00105 
00106 
00107 
00108 
00109 
00119 
00120 
00121 
00126 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00127 class AbstractNonlinearAssemblerSolverHybrid : public AbstractFeObjectAssembler<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM,true,true,NONLINEAR>
00128 {
00129 protected:
00131     AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* mpMesh;
00132 
00134     BoundaryConditionsContainer<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>* mpBoundaryConditions;
00135 
00137     AbstractNonlinearSolver* mpSolver;
00138 
00140     bool mWeAllocatedSolverMemory;
00141 
00143     bool mUseAnalyticalJacobian;
00144 
00145 
00146 
00156     void ApplyDirichletConditions(Vec currentGuess, Vec residual, Mat* pMat);
00157 
00165     void ComputeJacobianNumerically(const Vec currentGuess, Mat* pJacobian);
00166 
00167 public :
00168 
00178     void ComputeResidual(const Vec currentGuess, Vec residualVector);
00179 
00191     void ComputeJacobian(const Vec currentGuess, Mat* pJacobian);
00192 
00193 
00194 
00202     AbstractNonlinearAssemblerSolverHybrid(AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>* pMesh,
00203                                            BoundaryConditionsContainer<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>* pBoundaryConditions,
00204                                            unsigned numQuadPoints = 2);
00205 
00209     virtual ~AbstractNonlinearAssemblerSolverHybrid();
00210 
00219     virtual Vec Solve(Vec initialGuess, bool useAnalyticalJacobian=false);
00220 
00228     void SetNonlinearSolver(AbstractNonlinearSolver* pNonlinearSolver);
00229     
00230     
00244     bool VerifyJacobian(double tol);
00245 };
00246 
00247 
00248 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00249 AbstractNonlinearAssemblerSolverHybrid<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::AbstractNonlinearAssemblerSolverHybrid(
00250             AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>* pMesh,
00251             BoundaryConditionsContainer<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>* pBoundaryConditions,
00252             unsigned numQuadPoints)
00253     :  AbstractFeObjectAssembler<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM,true,true,NONLINEAR>(pMesh,numQuadPoints),
00254        mpMesh(pMesh),
00255        mpBoundaryConditions(pBoundaryConditions)
00256 {
00257     // if this is run with SPACE_DIM != ELEMENT_DIM the class has to be checked -
00258     // lots of places where should be using SPACE_DIM not ELEMENT_DIM??
00259     assert(SPACE_DIM==ELEMENT_DIM);
00260     assert(pMesh!=NULL);
00261     assert(pBoundaryConditions!=NULL);
00262 
00263    // mpAssembler = new SimpleNonlinearEllipticPdeAssembler<ELEMENT_DIM,SPACE_DIM>(mpMesh,pPde,numQuadPoints);
00264     mpSolver = new SimplePetscNonlinearSolver;
00265     mWeAllocatedSolverMemory = true;
00266 
00267     assert(mpMesh->GetNumNodes() == mpMesh->GetDistributedVectorFactory()->GetProblemSize());
00268 }
00269 
00270 
00271 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00272 AbstractNonlinearAssemblerSolverHybrid<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::~AbstractNonlinearAssemblerSolverHybrid()
00273 {
00274     if (mWeAllocatedSolverMemory)
00275     {
00276         delete mpSolver;
00277     }
00278 }
00279 
00280 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00281 void AbstractNonlinearAssemblerSolverHybrid<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::ApplyDirichletConditions(Vec currentGuess, Vec residual, Mat* pJacobian)
00282 {
00283     if (residual)
00284     {
00285         mpBoundaryConditions->ApplyDirichletToNonlinearResidual( currentGuess, 
00286                                                                  residual, 
00287                                                                  *(mpMesh->GetDistributedVectorFactory()));
00288     }
00289     if (pJacobian)
00290     {
00291         mpBoundaryConditions->ApplyDirichletToNonlinearJacobian(*pJacobian);
00292     }
00293 }
00294 
00295 
00296 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00297 void AbstractNonlinearAssemblerSolverHybrid<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::ComputeResidual(const Vec currentGuess, Vec residualVector)
00298 {
00299     this->SetVectorToAssemble(residualVector,true);
00300     this->SetCurrentSolution(currentGuess);
00301     this->SetApplyNeummanBoundaryConditionsToVector(mpBoundaryConditions);
00302     this->AssembleVector();
00303 
00304     PetscVecTools::Assemble(residualVector);
00305 
00306     ApplyDirichletConditions(currentGuess, residualVector, NULL);
00307 }
00308 
00309 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00310 void AbstractNonlinearAssemblerSolverHybrid<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::ComputeJacobian(const Vec currentGuess, Mat* pJacobian)
00311 {
00312     if (mUseAnalyticalJacobian)
00313     {
00314         this->SetMatrixToAssemble(*pJacobian);
00315         this->SetCurrentSolution(currentGuess);
00316         this->AssembleMatrix();
00317     
00318         PetscMatTools::AssembleIntermediate(*pJacobian);
00319         
00320         ApplyDirichletConditions(currentGuess, NULL, pJacobian);
00321 
00322         PetscMatTools::AssembleFinal(*pJacobian);
00323     }
00324     else
00325     {
00326         ComputeJacobianNumerically(currentGuess, pJacobian);
00327     }
00328 }
00329 
00330 
00331 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00332 void AbstractNonlinearAssemblerSolverHybrid<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::ComputeJacobianNumerically(const Vec currentGuess, Mat* pJacobian)
00333 {
00334     unsigned num_unknowns = PROBLEM_DIM * this->mpMesh->GetNumNodes();
00335 
00336     // Set up working vectors
00337     Vec residual;
00338     Vec perturbed_residual;
00339     Vec result;
00340     residual = PetscTools::CreateVec(num_unknowns);
00341     result = PetscTools::CreateVec(num_unknowns);
00342     perturbed_residual = PetscTools::CreateVec(num_unknowns);
00343     
00344     // Copy the currentGuess vector; we perturb the copy
00345     Vec current_guess_copy;
00346     PETSCEXCEPT( VecDuplicate(currentGuess, &current_guess_copy) );
00347     PETSCEXCEPT( VecCopy(currentGuess, current_guess_copy) );
00348 
00349     // Compute the current residual
00350     ComputeResidual(currentGuess, residual);
00351 
00352     // Amount to perturb each input element by
00353     double h = 1e-5;
00354 
00355     PetscInt ilo, ihi;
00356     VecGetOwnershipRange(current_guess_copy, &ilo, &ihi);
00357     unsigned lo = ilo;
00358     unsigned hi = ihi;
00359 
00360     // Iterate over entries in the input vector.
00361     for (unsigned global_index_outer = 0; global_index_outer < num_unknowns; global_index_outer++)
00362     {
00363         //Only perturb if we own it
00364         PetscVecTools::AddToElement(current_guess_copy, global_index_outer, h);
00365         ComputeResidual(current_guess_copy, perturbed_residual);
00366 
00367         // result = (perturbed_residual - residual) / h
00368         PetscVecTools::WAXPY(result, -1.0, residual, perturbed_residual);
00369         PetscVecTools::Scale(result, 1.0/h);
00370 
00371         double* p_result;
00372         PETSCEXCEPT( VecGetArray(result, &p_result) );
00373         for (unsigned global_index=lo; global_index < hi; global_index++)
00374         {
00375             unsigned local_index = global_index - lo;
00376             PetscMatTools::SetElement(*pJacobian, global_index, global_index_outer, p_result[local_index]);
00377         }
00378         PETSCEXCEPT( VecRestoreArray(result, &p_result) );
00379 
00380         PetscVecTools::AddToElement(current_guess_copy, global_index_outer, -h);
00381     }
00382 
00383     VecDestroy(residual);
00384     VecDestroy(perturbed_residual);
00385     VecDestroy(result);
00386     VecDestroy(current_guess_copy);
00387 
00388     PetscMatTools::AssembleFinal(*pJacobian);
00389 }
00390 
00391 
00392 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00393 Vec AbstractNonlinearAssemblerSolverHybrid<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::Solve(Vec initialGuess,
00394                                                                                        bool useAnalyticalJacobian)
00395 {
00396     assert(initialGuess!=NULL);
00397     mUseAnalyticalJacobian = useAnalyticalJacobian;
00398 
00399     PetscInt size_of_init_guess;
00400     VecGetSize(initialGuess, &size_of_init_guess);
00401     PetscInt problem_size = PROBLEM_DIM * this->mpMesh->GetNumNodes();
00402     if (size_of_init_guess != problem_size)
00403     {
00404         std::stringstream error_message;
00405         error_message << "Size of initial guess vector, " << size_of_init_guess
00406                       << ", does not match size of problem, " << problem_size;
00407         EXCEPTION(error_message.str());
00408     }
00409 
00410     // run the solver, telling it which global functions to call in order to assemble
00411     // the residual or jacobian
00412     return mpSolver->Solve(&AbstractNonlinearAssemblerSolverHybrid_ComputeResidual<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>,
00413                            &AbstractNonlinearAssemblerSolverHybrid_ComputeJacobian<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>,
00414                            initialGuess,
00415                            PROBLEM_DIM * this->mpMesh->CalculateMaximumNodeConnectivityPerProcess(), 
00416                            this);
00417 }
00418 
00419 
00420 
00421 
00422 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00423 void AbstractNonlinearAssemblerSolverHybrid<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::SetNonlinearSolver(AbstractNonlinearSolver* pNonlinearSolver)
00424 {
00425     if (mWeAllocatedSolverMemory)
00426     {
00427         delete mpSolver;
00428     }
00429     mpSolver = pNonlinearSolver;
00430     mWeAllocatedSolverMemory = false;
00431 }
00432 
00433 
00434 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00435 bool AbstractNonlinearAssemblerSolverHybrid<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::VerifyJacobian(double tol)
00436 {
00437     unsigned size = PROBLEM_DIM * this->mpMesh->GetNumNodes();
00438 
00439     Vec initial_guess = PetscTools::CreateAndSetVec(size, 0.0);
00440     
00441     Mat analytic_jacobian; 
00442     Mat numerical_jacobian;
00443 
00444     PetscTools::SetupMat(analytic_jacobian, size, size, size);
00445     PetscTools::SetupMat(numerical_jacobian, size, size, size);
00446 
00447     mUseAnalyticalJacobian = true;
00448     ComputeJacobian(initial_guess, &analytic_jacobian);
00449 
00450     mUseAnalyticalJacobian = false;
00451     ComputeJacobian(initial_guess, &numerical_jacobian);
00452     
00453     double minus_one = -1.0;
00454 #if (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 2) //PETSc 2.2
00455     // MatAYPX(*a, X, Y) does  Y = X + a*Y. 
00456     MatAYPX(&minus_one, analytic_jacobian, numerical_jacobian);
00457 #elif (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 3 && PETSC_VERSION_SUBMINOR == 1) //PETSc 2.3.1
00458     // MatAYPX( Y, a, X, structure) does Y = a*Y + X. 
00459     MatAYPX(numerical_jacobian, minus_one, analytic_jacobian); 
00460 #else
00461     // MatAYPX( Y, a, X, structure) does Y = a*Y + X. 
00462     MatAYPX(numerical_jacobian, minus_one, analytic_jacobian, DIFFERENT_NONZERO_PATTERN);
00463 #endif
00464     double norm;
00465     MatNorm(numerical_jacobian,NORM_INFINITY,&norm);
00466     
00467     MatDestroy(numerical_jacobian);
00468     MatDestroy(analytic_jacobian);
00469     VecDestroy(initial_guess);
00470 
00471     return (norm<tol);
00472 }
00473 
00474 
00475 
00476 
00486 
00487 
00488 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00489 PetscErrorCode AbstractNonlinearAssemblerSolverHybrid_ComputeResidual(SNES snes, Vec currentGuess,
00490                                                                       Vec residualVector, void* pContext)
00491 {
00492     // Extract the solver from the void*
00493     AbstractNonlinearAssemblerSolverHybrid<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>* p_solver =
00494         (AbstractNonlinearAssemblerSolverHybrid<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>*) pContext;
00495 
00496     p_solver->ComputeResidual(currentGuess, residualVector);
00497 
00498     return 0;
00499 }
00500 
00501 
00502 
00503 
00504 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00505 PetscErrorCode AbstractNonlinearAssemblerSolverHybrid_ComputeJacobian(SNES snes, Vec currentGuess,
00506                                                                       Mat* pGlobalJacobian, Mat* pPreconditioner,
00507                                                                       MatStructure* pMatStructure, void* pContext)
00508 {
00509     // Extract the solver from the void*
00510     AbstractNonlinearAssemblerSolverHybrid<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>* p_solver =
00511         (AbstractNonlinearAssemblerSolverHybrid<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>*) pContext;
00512 
00513     p_solver->ComputeJacobian(currentGuess, pGlobalJacobian);
00514 
00515     return 0;
00516 }
00517 
00518 
00519 #endif /*ABSTRACTNONLINEARASSEMBLERSOLVERHYBRID_HPP_*/

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