AbstractNonlinearAssembler.hpp

/*

Copyright (C) University of Oxford, 2005-2009

University of Oxford means the Chancellor, Masters and Scholars of the
University of Oxford, having an administrative office at Wellington
Square, Oxford OX1 2JD, UK.

This file is part of Chaste.

Chaste is free software: you can redistribute it and/or modify it
under the terms of the GNU Lesser General Public License as published
by the Free Software Foundation, either version 2.1 of the License, or
(at your option) any later version.

Chaste is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
License for more details. The offer of Chaste under the terms of the
License is subject to the License being interpreted in accordance with
English Law and subject to any action against the University of Oxford
being under the jurisdiction of the English Courts.

You should have received a copy of the GNU Lesser General Public License
along with Chaste. If not, see <http://www.gnu.org/licenses/>.

*/
#ifndef _ABSTRACTNONLINEARASSEMBLER_HPP_
#define _ABSTRACTNONLINEARASSEMBLER_HPP_

#include <vector>

#include "AbstractStaticAssembler.hpp"
#include "BoundaryConditionsContainer.hpp"
#include "AbstractNonlinearSolver.hpp"
#include "AbstractNonlinearSolver.hpp"
#include "SimplePetscNonlinearSolver.hpp"
#include "RandomNumberGenerator.hpp"
#include "PetscTools.hpp"



/*
 * Since we need to pass function pointers to the PETSc SNES routines, we can't
 * make these functions below methods. This is a pain, since it also means we
 * need to pass round a pointer to our assembler object as the void *pContext,
 * and cast it within the function to access data members.
 *
 * All the functions are defined as stubs which call methods on *pContext.
 *
 * Note: these are global functions, hence the need for long names to avoid
 * potential conflicting names later
 *
 * [The implementations are at the bottom of this file]
 */
template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM, class CONCRETE>
PetscErrorCode AbstractNonlinearAssembler_AssembleResidual(SNES snes,
        Vec currentGuess,
        Vec residualVector,
        void *pContext);

template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM, class CONCRETE>
PetscErrorCode AbstractNonlinearAssembler_AssembleJacobian(SNES snes,
        Vec currentGuess,
        Mat *pGlobalJacobian,
        Mat *pPreconditioner,
        MatStructure *pMatStructure,
        void *pContext);



template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM, class CONCRETE>
class AbstractNonlinearAssembler : public AbstractStaticAssembler<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM,true,CONCRETE>
{
    typedef AbstractStaticAssembler<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM,true,CONCRETE> BaseClassType;
private:
    Vec mInitialGuess;
protected:
    AbstractNonlinearSolver* mpSolver;
    bool mWeAllocatedSolverMemory;

    bool mUseAnalyticalJacobian;

    void ApplyDirichletConditions(Vec currentGuess, bool applyToMatrix)
    {
        Vec& residual = this->mpLinearSystem->rGetRhsVector();
        Mat& jacobian = this->mpLinearSystem->rGetLhsMatrix();
        assert((jacobian && applyToMatrix) || (!jacobian && !applyToMatrix));

        if (residual)
        {
            this->mpBoundaryConditions->ApplyDirichletToNonlinearResidual(currentGuess, residual);
        }
        if (jacobian)
        {
            this->mpBoundaryConditions->ApplyDirichletToNonlinearJacobian(jacobian);
        }
    }

public :
    PetscErrorCode AssembleResidual(const Vec currentGuess, Vec residualVector)
    {
        // call assemble system with the current guess and the residual vector
        // to be assembled
        this->PrepareForSolve();
        delete this->mpLinearSystem;
        this->mpLinearSystem = new LinearSystem(residualVector, NULL);
        this->AssembleSystem(true, false, currentGuess, 0.0);
        return 0;
    }


    PetscErrorCode AssembleJacobian(const Vec currentGuess, Mat *pGlobalJacobian)
    {
        if (mUseAnalyticalJacobian)
        {
            // call assemble system with the current guess and the jacobian to
            // be assembled
            delete this->mpLinearSystem;
            this->mpLinearSystem = new LinearSystem(NULL, *pGlobalJacobian);
            this->AssembleSystem(false, true, currentGuess, 0.0);
            return 0; // no error
        }
        else
        {
            return AssembleJacobianNumerically(currentGuess, pGlobalJacobian);
        }
    }

protected:
    PetscErrorCode AssembleJacobianNumerically(const Vec currentGuess, Mat *pJacobian)
    {
        unsigned num_nodes = PROBLEM_DIM*this->mpMesh->GetNumNodes();

        // Set up working vectors
        Vec residual;
        Vec perturbed_residual;
        Vec result;

        VecCreate(PETSC_COMM_WORLD, &residual);
        VecCreate(PETSC_COMM_WORLD, &result);
        VecCreate(PETSC_COMM_WORLD, &perturbed_residual);

        VecSetSizes(residual,          PETSC_DECIDE,num_nodes);
        VecSetSizes(result,            PETSC_DECIDE,num_nodes);
        VecSetSizes(perturbed_residual,PETSC_DECIDE,num_nodes);

        VecSetFromOptions(residual);
        VecSetFromOptions(result);
        VecSetFromOptions(perturbed_residual);

        // Copy the currentGuess vector; we perturb the copy
        Vec current_guess_copy;
        PETSCEXCEPT( VecDuplicate(currentGuess, &current_guess_copy) );
        PETSCEXCEPT( VecCopy(currentGuess, current_guess_copy) );

        // Compute the current residual
        AssembleResidual(currentGuess, residual);

        // Amount to perturb each input element by
        double h = 0.00001;
        PetscScalar subtract = -1;
        PetscScalar one_over_h = 1.0/h;

        PetscInt ilo, ihi;
        VecGetOwnershipRange(current_guess_copy, &ilo, &ihi);
        unsigned lo=ilo;
        unsigned hi=ihi;

        // Iterate over entries in the input vector.
        for (unsigned global_index_outer = 0; global_index_outer < num_nodes; global_index_outer++)
        {
            //Only perturb if we own it
            if (lo<=global_index_outer && global_index_outer<hi)
            {
                PETSCEXCEPT( VecSetValue(current_guess_copy, global_index_outer,h, ADD_VALUES) );
            }
            AssembleResidual(current_guess_copy, perturbed_residual);

            // result = (perturbed_residual - residual) / h
#if (PETSC_VERSION_MINOR == 2) //Old API
            PETSCEXCEPT( VecWAXPY(&subtract, residual, perturbed_residual, result) );
            PETSCEXCEPT( VecScale(&one_over_h, result) );
#else
            PETSCEXCEPT( VecWAXPY(result, subtract, residual, perturbed_residual) );
            PETSCEXCEPT( VecScale(result, one_over_h) );
#endif

            double *p_result;
            PETSCEXCEPT( VecGetArray(result, &p_result) );
            for (unsigned global_index=lo; global_index < hi; global_index++)
            {
                unsigned local_index = global_index - lo;
                PETSCEXCEPT( MatSetValue(*pJacobian, global_index, global_index_outer,
                                         p_result[local_index], INSERT_VALUES) );
            }
            PETSCEXCEPT( VecRestoreArray(result, &p_result) );

            if (lo<=global_index_outer && global_index_outer<hi)
            {
                PETSCEXCEPT( VecSetValue(current_guess_copy, global_index_outer, -h, ADD_VALUES) );
            }
        }

        VecDestroy(residual);
        VecDestroy(perturbed_residual);
        VecDestroy(result);
        VecDestroy(current_guess_copy);

        MatAssemblyBegin(*pJacobian, MAT_FINAL_ASSEMBLY);
        MatAssemblyEnd(*pJacobian, MAT_FINAL_ASSEMBLY);

        return 0; // No error
    }

    virtual void AssembleSystem(bool assembleVector, bool assembleMatrix,
                                Vec currentGuess=NULL, double currentTime=0.0)
    {
        // If the problem is nonlinear the currentGuess MUST be specifed
        assert( currentGuess );
        BaseClassType::AssembleSystem(assembleVector, assembleMatrix, currentGuess, currentTime);
    }

    bool ProblemIsNonlinear()
    {
        return true;
    }

    void InitialiseForSolve(Vec initialGuess)
    {
        // Check size of initial guess is correct
        PetscInt size_of_init_guess;
        VecGetSize(initialGuess, &size_of_init_guess);
        PetscInt problem_size = PROBLEM_DIM * this->mpMesh->GetNumNodes();
        if (size_of_init_guess != problem_size)
        {
            std::stringstream error_message;
            error_message << "Size of initial guess vector, " << size_of_init_guess
                          << ", does not match size of problem, " << problem_size;
            EXCEPTION(error_message.str());
        }
    }

    Vec StaticSolve(Vec currentSolutionOrGuess=NULL,
                    double currentTime=0.0,
                    bool assembleMatrix=true)
    {
        assert(this->mpBoundaryConditions!=NULL);
        assert(currentSolutionOrGuess!=NULL);
        assert(assembleMatrix); 

        // run the solver, telling it which global functions to call in order to assemble
        // the residual or jacobian
        Vec answer = this->mpSolver->Solve(&AbstractNonlinearAssembler_AssembleResidual<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM, CONCRETE>,
                                           &AbstractNonlinearAssembler_AssembleJacobian<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM, CONCRETE>,
                                           currentSolutionOrGuess,
                                           this);
        return answer;
    }

public:
    AbstractNonlinearAssembler(unsigned numQuadPoints = 2) :
            BaseClassType(numQuadPoints)
    {
        mpSolver = new SimplePetscNonlinearSolver;
        mWeAllocatedSolverMemory = true;

        mUseAnalyticalJacobian = false;

        mInitialGuess = NULL;
    }

    ~AbstractNonlinearAssembler()
    {
        if (mWeAllocatedSolverMemory)
        {
            delete mpSolver;
        }
    }

    void SetUseAnalyticalJacobian(bool useAnalyticalJacobian)
    {
        mUseAnalyticalJacobian = useAnalyticalJacobian;
    }

    virtual Vec Solve(Vec initialGuess, bool useAnalyticalJacobian=false)
    {
        assert(initialGuess!=NULL);
        SetUseAnalyticalJacobian(useAnalyticalJacobian);

        this->PrepareForSolve();
        InitialiseForSolve(initialGuess);

        return StaticSolve(initialGuess);
    }

    void SetNonlinearSolver(AbstractNonlinearSolver* pNonlinearSolver)
    {
        if (mWeAllocatedSolverMemory)
        {
            delete mpSolver;
        }
        mpSolver = pNonlinearSolver;
        mWeAllocatedSolverMemory = false;
    }

    Vec CreateConstantInitialGuess(double value)
    {
        assert(this->mpMesh!=NULL);
        unsigned size = PROBLEM_DIM * this->mpMesh->GetNumNodes();
        return PetscTools::CreateVec(size, value);
    }

    bool VerifyJacobian(double tol=1e-4, bool print=false)
    {
        unsigned size = PROBLEM_DIM * this->mpMesh->GetNumNodes();

        Vec initial_guess;
        VecCreate(PETSC_COMM_WORLD, &initial_guess);
        VecSetSizes(initial_guess, PETSC_DECIDE, size);
        VecSetFromOptions(initial_guess);


        for (unsigned i=0; i<size; i++)
        {
            VecSetValue(initial_guess, i, 0.0, INSERT_VALUES);
        }

        VecAssemblyBegin(initial_guess);
        VecAssemblyEnd(initial_guess);


        Mat analytic_jacobian; //Jacobian Matrix
        Mat numerical_jacobian; //Jacobian Matrix

        PetscTools::SetupMat(analytic_jacobian, size, size);
        PetscTools::SetupMat(numerical_jacobian, size, size);

        mUseAnalyticalJacobian = true;
        AssembleJacobian(initial_guess, &analytic_jacobian);

        mUseAnalyticalJacobian = false;
        AssembleJacobian(initial_guess, &numerical_jacobian);

        bool all_less_than_tol = true;

        if (print)
        {
            std::cout << "Difference between numerical and analyical Jacobians:\n\n";
        }
        for (unsigned i=0; i<size; i++)
        {
            for (unsigned j=0; j<size; j++)
            {
                double val_a[1];
                double val_n[1];
                PetscInt row[1];
                PetscInt col[1];
                row[0] = i;
                col[0] = j;

                MatGetValues(numerical_jacobian,1,row,1,col,val_n);
                MatGetValues(analytic_jacobian,1,row,1,col,val_a);

                if (print)
                {
                    std::cout << val_n[0] - val_a[0]<< " ";
                }

                if (fabs(val_n[0]-val_a[0]) > tol)
                {
#define COVERAGE_IGNORE // would have to write a bad concrete assembler class just to cover this line
                    all_less_than_tol = false;
#undef COVERAGE_IGNORE
                }
            }
            if (print)
            {
                std::cout << "\n" <<std::flush;
            }
        }
        std::cout << std::flush;

        MatDestroy(numerical_jacobian);
        MatDestroy(analytic_jacobian);
        VecDestroy(initial_guess);

        return all_less_than_tol;
    }

}; //end of class AbstractNonlinearAssembler






//========================================================================================//
//========================================================================================//
//                                                                                        //
//                Global functions called by the Petsc nonlinear solver                   //
//                                                                                        //
//========================================================================================//
//========================================================================================//



template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM, class CONCRETE>
PetscErrorCode AbstractNonlinearAssembler_AssembleResidual(SNES snes, Vec currentGuess,
        Vec residualVector, void *pContext)
{
    // Extract an assembler from the void*
    AbstractNonlinearAssembler<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM, CONCRETE> *pAssembler =
        (AbstractNonlinearAssembler<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM, CONCRETE>*) pContext;

    PetscErrorCode ierr = pAssembler->AssembleResidual(currentGuess, residualVector);

    //double two_norm;
    //VecNorm(residualVector, NORM_2, &two_norm);
    //std::cout << "||residual|| = " << two_norm << "\n";

    return ierr;
}



template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM, class CONCRETE>
PetscErrorCode AbstractNonlinearAssembler_AssembleJacobian(SNES snes, Vec currentGuess,
        Mat *pGlobalJacobian, Mat *pPreconditioner,
        MatStructure *pMatStructure, void *pContext)
{
    //std::cout << "begin assemble jacobian\n";

    // Extract an assembler from the void*
    AbstractNonlinearAssembler<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM, CONCRETE> *pAssembler =
        (AbstractNonlinearAssembler<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM, CONCRETE>*) pContext;

    PetscErrorCode ierr = pAssembler->AssembleJacobian(currentGuess, pGlobalJacobian);
    //std::cout << "end assemble jacobian\n";

    return ierr;
}



#endif //_ABSTRACTNONLINEARASSEMBLER_HPP_

---