Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
AbstractNonlinearAssemblerSolverHybrid.hpp
1/*
2
3Copyright (c) 2005-2024, University of Oxford.
4All rights reserved.
5
6University of Oxford means the Chancellor, Masters and Scholars of the
7University of Oxford, having an administrative office at Wellington
8Square, Oxford OX1 2JD, UK.
9
10This file is part of Chaste.
11
12Redistribution and use in source and binary forms, with or without
13modification, are permitted provided that the following conditions are met:
14 * Redistributions of source code must retain the above copyright notice,
15 this list of conditions and the following disclaimer.
16 * Redistributions in binary form must reproduce the above copyright notice,
17 this list of conditions and the following disclaimer in the documentation
18 and/or other materials provided with the distribution.
19 * Neither the name of the University of Oxford nor the names of its
20 contributors may be used to endorse or promote products derived from this
21 software without specific prior written permission.
22
23THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33
34*/
35
36#ifndef ABSTRACTNONLINEARASSEMBLERSOLVERHYBRID_HPP_
37#define ABSTRACTNONLINEARASSEMBLERSOLVERHYBRID_HPP_
38
39#include "BoundaryConditionsContainer.hpp"
40#include "AbstractNonlinearEllipticPde.hpp"
41#include "AbstractNonlinearSolver.hpp"
42#include "PetscTools.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"
49
50/*
51 * Definitions of GLOBAL functions used by PETSc nonlinear solver
52 * (implementations are at the bottom of this file).
53 */
54
69template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
70PetscErrorCode AbstractNonlinearAssemblerSolverHybrid_ComputeResidual(SNES snes,
71 Vec currentGuess,
72 Vec residualVector,
73 void* pContext);
74#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>=5)
93 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
94
95 PetscErrorCode AbstractNonlinearAssemblerSolverHybrid_ComputeJacobian(SNES snes,
96 Vec currentGuess,
97 Mat globalJacobian,
98 Mat preconditioner,
99 void* pContext);
100#else
117 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
118 PetscErrorCode AbstractNonlinearAssemblerSolverHybrid_ComputeJacobian(SNES snes,
119 Vec currentGuess,
120 Mat* pGlobalJacobian,
121 Mat* pPreconditioner,
122 MatStructure* pMatStructure,
123 void* pContext);
124#endif
125
131template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
132class AbstractNonlinearAssemblerSolverHybrid : public AbstractFeVolumeIntegralAssembler<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM,true,true,NONLINEAR>
133{
134protected:
135
138
141
144
147
150
153
163 void ApplyDirichletConditions(Vec currentGuess, Vec residual, Mat* pMat);
164
172 void ComputeJacobianNumerically(const Vec currentGuess, Mat* pJacobian);
173
174public:
175
185 void ComputeResidual(const Vec currentGuess, Vec residualVector);
186
198 void ComputeJacobian(const Vec currentGuess, Mat* pJacobian);
199
208
213
223 virtual Vec Solve(Vec initialGuess, bool useAnalyticalJacobian=false);
224
233
247 bool VerifyJacobian(double tol);
248};
249
250template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
254 : AbstractFeVolumeIntegralAssembler<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM,true,true,NONLINEAR>(pMesh),
255 mpMesh(pMesh),
256 mpBoundaryConditions(pBoundaryConditions)
257{
258 /*
259 * Note: if this is run with SPACE_DIM != ELEMENT_DIM the class has to be checked:
260 * there may be lots of places where we should be using SPACE_DIM not ELEMENT_DIM.
261 */
262 assert(SPACE_DIM==ELEMENT_DIM);
263 assert(pMesh!=nullptr);
264 assert(pBoundaryConditions!=nullptr);
265
268
269 assert(mpMesh->GetNumNodes() == mpMesh->GetDistributedVectorFactory()->GetProblemSize());
270
272}
273
274template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
276{
277 if (mWeAllocatedSolverMemory)
278 {
279 delete mpSolver;
280 }
281 delete mpNeumannSurfaceTermsAssembler;
282}
283
284template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
286{
287 if (residual)
288 {
289 mpBoundaryConditions->ApplyDirichletToNonlinearResidual(currentGuess,
290 residual,
291 *(mpMesh->GetDistributedVectorFactory()));
292 }
293 if (pJacobian)
294 {
295 mpBoundaryConditions->ApplyDirichletToNonlinearJacobian(*pJacobian);
296 }
297}
298
299template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
301{
302 // Add the volume integral part of the residual. This will call ComputeVectorTerm, which needs to be
303 // implemented in the concrete class.
304 this->SetVectorToAssemble(residualVector,true);
305 this->SetCurrentSolution(currentGuess);
306 this->AssembleVector();
307
308 // Add the surface integral contribution - note the negative sign (scale factor).
309 // This contribution is, for example for a 1 unknown problem
310 // -integral(g\phi_id dS), where g is the specfied neumann boundary condition
311 mpNeumannSurfaceTermsAssembler->SetVectorToAssemble(residualVector, false);
312 mpNeumannSurfaceTermsAssembler->SetScaleFactor(-1.0);
313 mpNeumannSurfaceTermsAssembler->Assemble();
314
315 PetscVecTools::Finalise(residualVector);
316
317 ApplyDirichletConditions(currentGuess, residualVector, nullptr);
318}
319
320template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
322{
323 if (mUseAnalyticalJacobian)
324 {
325 this->SetMatrixToAssemble(*pJacobian);
326 this->SetCurrentSolution(currentGuess);
327 this->AssembleMatrix();
328
330
331 ApplyDirichletConditions(currentGuess, nullptr, pJacobian);
332
333 PetscMatTools::Finalise(*pJacobian);
334 }
335 else
336 {
337 ComputeJacobianNumerically(currentGuess, pJacobian);
338 }
339}
340
341template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
343{
344 unsigned num_unknowns = PROBLEM_DIM * this->mpMesh->GetNumNodes();
345
346 // Set up working vectors
347 Vec residual;
348 Vec perturbed_residual;
349 Vec result;
350 residual = PetscTools::CreateVec(num_unknowns);
351 result = PetscTools::CreateVec(num_unknowns);
352 perturbed_residual = PetscTools::CreateVec(num_unknowns);
353
354 // Copy the currentGuess vector; we perturb the copy
355 Vec current_guess_copy;
356 PETSCEXCEPT( VecDuplicate(currentGuess, &current_guess_copy) );
357 PETSCEXCEPT( VecCopy(currentGuess, current_guess_copy) );
358
359 // Compute the current residual
360 ComputeResidual(currentGuess, residual);
361
362 // Amount to perturb each input element by
363 double h = 1e-5;
364 double near_hsquared = 1e-9;
365
366 PetscInt ilo, ihi;
367 VecGetOwnershipRange(current_guess_copy, &ilo, &ihi);
368 unsigned lo = ilo;
369 unsigned hi = ihi;
370
371 // Iterate over entries in the input vector
372 for (unsigned global_index_outer = 0; global_index_outer < num_unknowns; global_index_outer++)
373 {
374 // Only perturb if we own it
375 PetscVecTools::AddToElement(current_guess_copy, global_index_outer, h);
376 ComputeResidual(current_guess_copy, perturbed_residual);
377
378 // result = (perturbed_residual - residual) / h
379 PetscVecTools::WAXPY(result, -1.0, residual, perturbed_residual);
380 PetscVecTools::Scale(result, 1.0/h);
381
382 double* p_result;
385 PETSCEXCEPT( VecGetArray(result, &p_result) );
386 for (unsigned global_index=lo; global_index < hi; global_index++)
387 {
388 double result_entry = p_result[ global_index - lo];
389 if (!CompareDoubles::IsNearZero(result_entry, near_hsquared))
390 {
391 PetscMatTools::SetElement(*pJacobian, global_index, global_index_outer, result_entry);
392 }
393 }
394 PETSCEXCEPT( VecRestoreArray(result, &p_result) );
395
396 PetscVecTools::AddToElement(current_guess_copy, global_index_outer, -h);
397 }
398
399 PetscTools::Destroy(residual);
400 PetscTools::Destroy(perturbed_residual);
401 PetscTools::Destroy(result);
402 PetscTools::Destroy(current_guess_copy);
403
404 PetscMatTools::Finalise(*pJacobian);
405}
406
407template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
409 bool useAnalyticalJacobian)
410{
411 assert(initialGuess != nullptr);
412 mUseAnalyticalJacobian = useAnalyticalJacobian;
413
414 PetscInt size_of_init_guess;
415 VecGetSize(initialGuess, &size_of_init_guess);
416 PetscInt problem_size = PROBLEM_DIM * this->mpMesh->GetNumNodes();
417 if (size_of_init_guess != problem_size)
418 {
419 EXCEPTION("Size of initial guess vector, " << size_of_init_guess
420 << ", does not match size of problem, " << problem_size);
421 }
422
423 // Run the solver, telling it which global functions to call in order to assemble the residual or jacobian
424 return mpSolver->Solve(&AbstractNonlinearAssemblerSolverHybrid_ComputeResidual<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>,
425 &AbstractNonlinearAssemblerSolverHybrid_ComputeJacobian<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>,
426 initialGuess,
427 PROBLEM_DIM * this->mpMesh->CalculateMaximumNodeConnectivityPerProcess(),
428 this);
429}
430
431template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
433{
434 if (mWeAllocatedSolverMemory)
435 {
436 delete mpSolver;
437 }
438 mpSolver = pNonlinearSolver;
439 mWeAllocatedSolverMemory = false;
440}
441
442template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
444{
445 unsigned size = PROBLEM_DIM * this->mpMesh->GetNumNodes();
446
447 Vec initial_guess = PetscTools::CreateAndSetVec(size, 0.0);
448
449 Mat analytic_jacobian;
450 Mat numerical_jacobian;
451
452 PetscTools::SetupMat(analytic_jacobian, size, size, size);
453 PetscTools::SetupMat(numerical_jacobian, size, size, size);
454
455 mUseAnalyticalJacobian = true;
456 ComputeJacobian(initial_guess, &analytic_jacobian);
457
458 mUseAnalyticalJacobian = false;
459 ComputeJacobian(initial_guess, &numerical_jacobian);
460
461 bool ok = PetscMatTools::CheckEquality(numerical_jacobian, analytic_jacobian, tol);
462 PetscTools::Destroy(numerical_jacobian);
463 PetscTools::Destroy(analytic_jacobian);
464 PetscTools::Destroy(initial_guess);
465
466 return ok;
467}
468
469// Implementations of GLOBAL functions called by PETSc nonlinear solver
470
471template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
472PetscErrorCode AbstractNonlinearAssemblerSolverHybrid_ComputeResidual(SNES snes,
473 Vec currentGuess,
474 Vec residualVector,
475 void* pContext)
476{
477 // Extract the solver from the void*
480
481 p_solver->ComputeResidual(currentGuess, residualVector);
482
483 return 0;
484}
485
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,
489 Vec currentGuess,
490 Mat globalJacobian,
491 Mat preconditioner,
492 void* pContext)
493{
494 // Extract the solver from the void*
497
498 p_solver->ComputeJacobian(currentGuess, &globalJacobian);
499
500 return 0;
501}
502#else
503PetscErrorCode AbstractNonlinearAssemblerSolverHybrid_ComputeJacobian(SNES snes,
504 Vec currentGuess,
505 Mat* pGlobalJacobian,
506 Mat* pPreconditioner,
507 MatStructure* pMatStructure,
508 void* pContext)
509{
510 // Extract the solver from the void*
513
514 p_solver->ComputeJacobian(currentGuess, pGlobalJacobian);
515
516 return 0;
517}
518#endif
519
520#endif /*ABSTRACTNONLINEARASSEMBLERSOLVERHYBRID_HPP_*/
#define EXCEPTION(message)
void ComputeJacobianNumerically(const Vec currentGuess, Mat *pJacobian)
void ComputeResidual(const Vec currentGuess, Vec residualVector)
void ComputeJacobian(const Vec currentGuess, Mat *pJacobian)
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
AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > * mpMesh
void ApplyDirichletConditions(Vec currentGuess, Vec residual, Mat *pMat)
NaturalNeumannSurfaceTermAssembler< ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM > * mpNeumannSurfaceTermsAssembler
void SetNonlinearSolver(AbstractNonlinearSolver *pNonlinearSolver)
static bool IsNearZero(double number, double tolerance)
static void SwitchWriteMode(Mat matrix)
static void Finalise(Mat matrix)
static void SetElement(Mat matrix, PetscInt row, PetscInt col, double value)
static bool CheckEquality(const Mat mat1, const Mat mat2, double tol=1e-10)
static Vec CreateAndSetVec(int size, double value)
static void Destroy(Vec &rVec)
static Vec CreateVec(int size, int localSize=PETSC_DECIDE, bool ignoreOffProcEntries=true)
static void SetupMat(Mat &rMat, int numRows, int numColumns, unsigned rowPreallocation, int numLocalRows=PETSC_DECIDE, int numLocalColumns=PETSC_DECIDE, bool ignoreOffProcEntries=true, bool newAllocationError=true)
static void WAXPY(Vec w, double a, Vec x, Vec y)
static void Scale(Vec vector, double scaleFactor)
static void AddToElement(Vec vector, PetscInt row, double value)
static void Finalise(Vec vector)