36#ifndef ABSTRACTNONLINEARASSEMBLERSOLVERHYBRID_HPP_
37#define ABSTRACTNONLINEARASSEMBLERSOLVERHYBRID_HPP_
39#include "BoundaryConditionsContainer.hpp"
40#include "AbstractNonlinearEllipticPde.hpp"
41#include "AbstractNonlinearSolver.hpp"
43#include "PetscMatTools.hpp"
44#include "PetscVecTools.hpp"
45#include "AbstractFeVolumeIntegralAssembler.hpp"
46#include "LinearSystem.hpp"
47#include "SimplePetscNonlinearSolver.hpp"
48#include "NaturalNeumannSurfaceTermAssembler.hpp"
69template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
70PetscErrorCode AbstractNonlinearAssemblerSolverHybrid_ComputeResidual(SNES snes,
74#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>=5)
93 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
95 PetscErrorCode AbstractNonlinearAssemblerSolverHybrid_ComputeJacobian(SNES snes,
117 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
118 PetscErrorCode AbstractNonlinearAssemblerSolverHybrid_ComputeJacobian(SNES snes,
120 Mat* pGlobalJacobian,
121 Mat* pPreconditioner,
122 MatStructure* pMatStructure,
131template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
223 virtual Vec Solve(
Vec initialGuess,
bool useAnalyticalJacobian=
false);
250template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
256 mpBoundaryConditions(pBoundaryConditions)
262 assert(SPACE_DIM==ELEMENT_DIM);
263 assert(pMesh!=
nullptr);
264 assert(pBoundaryConditions!=
nullptr);
269 assert(
mpMesh->GetNumNodes() ==
mpMesh->GetDistributedVectorFactory()->GetProblemSize());
274template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
277 if (mWeAllocatedSolverMemory)
281 delete mpNeumannSurfaceTermsAssembler;
284template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
289 mpBoundaryConditions->ApplyDirichletToNonlinearResidual(currentGuess,
291 *(mpMesh->GetDistributedVectorFactory()));
295 mpBoundaryConditions->ApplyDirichletToNonlinearJacobian(*pJacobian);
299template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
304 this->SetVectorToAssemble(residualVector,
true);
305 this->SetCurrentSolution(currentGuess);
306 this->AssembleVector();
311 mpNeumannSurfaceTermsAssembler->SetVectorToAssemble(residualVector,
false);
312 mpNeumannSurfaceTermsAssembler->SetScaleFactor(-1.0);
313 mpNeumannSurfaceTermsAssembler->Assemble();
317 ApplyDirichletConditions(currentGuess, residualVector,
nullptr);
320template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
323 if (mUseAnalyticalJacobian)
325 this->SetMatrixToAssemble(*pJacobian);
326 this->SetCurrentSolution(currentGuess);
327 this->AssembleMatrix();
331 ApplyDirichletConditions(currentGuess,
nullptr, pJacobian);
337 ComputeJacobianNumerically(currentGuess, pJacobian);
341template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
344 unsigned num_unknowns = PROBLEM_DIM * this->mpMesh->GetNumNodes();
348 Vec perturbed_residual;
355 Vec current_guess_copy;
356 PETSCEXCEPT( VecDuplicate(currentGuess, ¤t_guess_copy) );
357 PETSCEXCEPT( VecCopy(currentGuess, current_guess_copy) );
360 ComputeResidual(currentGuess, residual);
364 double near_hsquared = 1e-9;
367 VecGetOwnershipRange(current_guess_copy, &ilo, &ihi);
372 for (
unsigned global_index_outer = 0; global_index_outer < num_unknowns; global_index_outer++)
376 ComputeResidual(current_guess_copy, perturbed_residual);
385 PETSCEXCEPT( VecGetArray(result, &p_result) );
386 for (
unsigned global_index=lo; global_index < hi; global_index++)
388 double result_entry = p_result[ global_index - lo];
394 PETSCEXCEPT( VecRestoreArray(result, &p_result) );
407template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
409 bool useAnalyticalJacobian)
411 assert(initialGuess !=
nullptr);
412 mUseAnalyticalJacobian = useAnalyticalJacobian;
415 VecGetSize(initialGuess, &size_of_init_guess);
416 PetscInt problem_size = PROBLEM_DIM * this->mpMesh->GetNumNodes();
417 if (size_of_init_guess != problem_size)
419 EXCEPTION(
"Size of initial guess vector, " << size_of_init_guess
420 <<
", does not match size of problem, " << problem_size);
424 return mpSolver->Solve(&AbstractNonlinearAssemblerSolverHybrid_ComputeResidual<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>,
425 &AbstractNonlinearAssemblerSolverHybrid_ComputeJacobian<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>,
427 PROBLEM_DIM * this->mpMesh->CalculateMaximumNodeConnectivityPerProcess(),
431template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
434 if (mWeAllocatedSolverMemory)
438 mpSolver = pNonlinearSolver;
439 mWeAllocatedSolverMemory =
false;
442template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
445 unsigned size = PROBLEM_DIM * this->mpMesh->GetNumNodes();
449 Mat analytic_jacobian;
450 Mat numerical_jacobian;
455 mUseAnalyticalJacobian =
true;
456 ComputeJacobian(initial_guess, &analytic_jacobian);
458 mUseAnalyticalJacobian =
false;
459 ComputeJacobian(initial_guess, &numerical_jacobian);
471template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
472PetscErrorCode AbstractNonlinearAssemblerSolverHybrid_ComputeResidual(SNES snes,
486template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
487#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>=5)
488PetscErrorCode AbstractNonlinearAssemblerSolverHybrid_ComputeJacobian(SNES snes,
503PetscErrorCode AbstractNonlinearAssemblerSolverHybrid_ComputeJacobian(SNES snes,
505 Mat* pGlobalJacobian,
506 Mat* pPreconditioner,
507 MatStructure* pMatStructure,
#define EXCEPTION(message)
void ComputeJacobianNumerically(const Vec currentGuess, Mat *pJacobian)
void ComputeResidual(const Vec currentGuess, Vec residualVector)
void ComputeJacobian(const Vec currentGuess, Mat *pJacobian)
AbstractNonlinearSolver * mpSolver
bool mUseAnalyticalJacobian
virtual Vec Solve(Vec initialGuess, bool useAnalyticalJacobian=false)
AbstractNonlinearAssemblerSolverHybrid(AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > *pMesh, BoundaryConditionsContainer< ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM > *pBoundaryConditions)
BoundaryConditionsContainer< ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM > * mpBoundaryConditions
bool VerifyJacobian(double tol)
AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > * mpMesh
void ApplyDirichletConditions(Vec currentGuess, Vec residual, Mat *pMat)
bool mWeAllocatedSolverMemory
NaturalNeumannSurfaceTermAssembler< ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM > * mpNeumannSurfaceTermsAssembler
void SetNonlinearSolver(AbstractNonlinearSolver *pNonlinearSolver)
virtual ~AbstractNonlinearAssemblerSolverHybrid()
static bool IsNearZero(double number, double tolerance)