Skip to content

Commit

Permalink
apply bx norm factor for eigen kernel tagging idaholab#28984
Browse files Browse the repository at this point in the history
  • Loading branch information
YaqiWang committed Nov 1, 2024
1 parent c80f515 commit 7948ad5
Show file tree
Hide file tree
Showing 8 changed files with 108 additions and 48 deletions.
6 changes: 6 additions & 0 deletions framework/include/interfaces/TaggingInterface.h
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,7 @@ class MooseObject;
class SubProblem;
class Assembly;
class ReferenceResidualProblem;
class EigenProblem;

template <typename T>
InputParameters validParams();
Expand Down Expand Up @@ -395,6 +396,11 @@ class TaggingInterface
/// this data member to avoid constant dynamic heap allocations
std::vector<Real> _absolute_residuals;

/// Pointer to the eigen problem for an eigen object
const EigenProblem * _eigen_problem = nullptr;
/// Pointer to the postprocessor for the Bx norm
const PostprocessorValue * _bx_norm = nullptr;

friend void NonlinearSystemBase::constraintJacobians(bool);
};

Expand Down
13 changes: 12 additions & 1 deletion framework/include/problems/EigenProblem.h
Original file line number Diff line number Diff line change
Expand Up @@ -49,6 +49,7 @@ class EigenProblem : public FEProblemBase

NonlinearEigenSystem & getNonlinearEigenSystem(const unsigned int nl_sys_num);
NonlinearEigenSystem & getCurrentNonlinearEigenSystem();
const NonlinearEigenSystem & getCurrentNonlinearEigenSystem() const;

virtual void checkProblemIntegrity() override;

Expand Down Expand Up @@ -227,7 +228,7 @@ class EigenProblem : public FEProblemBase
/**
* Form the Bx norm
*/
Real formNorm();
const Real & formNorm() const;

/**
* Whether a Bx norm postprocessor has been provided
Expand All @@ -239,6 +240,8 @@ class EigenProblem : public FEProblemBase
*/
void setBxNorm(const PostprocessorName & bx_norm) { _bx_norm_name = bx_norm; }

const PostprocessorName & bxNormName() const { return _bx_norm_name.value(); }

protected:
unsigned int _n_eigen_pairs_required;
bool _generalized_eigenvalue_problem;
Expand Down Expand Up @@ -291,6 +294,8 @@ class EigenProblem : public FEProblemBase
/// default L2 norm of Bx will be used as the Bx norm
std::optional<PostprocessorName> _bx_norm_name;

/// Pointer to the postprocessor providing the Bx norm
const PostprocessorValue * _bx_norm = nullptr;
#endif

using FEProblemBase::_nl;
Expand All @@ -312,4 +317,10 @@ EigenProblem::getCurrentNonlinearEigenSystem()
return *_nl_eigen;
}

inline const NonlinearEigenSystem &
EigenProblem::getCurrentNonlinearEigenSystem() const
{
return *_nl_eigen;
}

#endif
29 changes: 29 additions & 0 deletions framework/src/interfaces/TaggingInterface.C
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,8 @@
#include "FEProblem.h"
#include "Assembly.h"
#include "ReferenceResidualProblem.h"
#include "EigenProblem.h"
#include "NonlinearEigenSystem.h"

#include "libmesh/dense_vector.h"

Expand Down Expand Up @@ -147,6 +149,19 @@ TaggingInterface::TaggingInterface(const MooseObject * moose_object)
_non_ref_vector_tags = _vector_tags;
_non_ref_abs_vector_tags = _abs_vector_tags;
}

// we need to scale the eigen kernel residual when tagging to a vector different than Bx vector
if (const auto * const eigen_problem = dynamic_cast<const EigenProblem *>(fe_problem))
{
const auto & eigen_system = eigen_problem->getCurrentNonlinearEigenSystem();
if (_vector_tags.count(eigen_system.eigenVectorTag()))
{
_eigen_problem = eigen_problem;

if (_eigen_problem->bxNormProvided())
_bx_norm = &_eigen_problem->getPostprocessorValueByName(eigen_problem->bxNormName());
}
}
}

void
Expand Down Expand Up @@ -359,8 +374,15 @@ TaggingInterface::prepareMatrixTagLower(Assembly & assembly,
void
TaggingInterface::accumulateTaggedLocalResidual()
{
if (_eigen_problem)
{
Real fac = (_eigen_problem->negativeSignEigenKernel() ? 1 : -1) * (_bx_norm ? *_bx_norm : 1);
_local_re.scale(1 / fac);
}

for (auto & re : _re_blocks)
*re += _local_re;

for (auto & absre : _absre_blocks)
for (const auto i : index_range(_local_re))
(*absre)(i) += std::abs(_local_re(i));
Expand All @@ -369,8 +391,15 @@ TaggingInterface::accumulateTaggedLocalResidual()
void
TaggingInterface::assignTaggedLocalResidual()
{
if (_eigen_problem)
{
Real fac = (_eigen_problem->negativeSignEigenKernel() ? 1 : -1) * (_bx_norm ? *_bx_norm : 1);
_local_re.scale(1 / fac);
}

for (auto & re : _re_blocks)
*re = _local_re;

for (auto & absre : _absre_blocks)
for (const auto i : index_range(_local_re))
(*absre)(i) = std::abs(_local_re(i));
Expand Down
25 changes: 18 additions & 7 deletions framework/src/problems/EigenProblem.C
Original file line number Diff line number Diff line change
Expand Up @@ -312,6 +312,9 @@ EigenProblem::computeResidualTag(const NumericVector<Number> & soln,
setCurrentNonlinearSystem(_nl_eigen->number());
computeResidualTags(_fe_vector_tags);

if (bxNormProvided() && tag == _nl_eigen->eigenVectorTag())
residual *= (-*_bx_norm);

_nl_eigen->disassociateVectorFromTag(residual, tag);
}

Expand Down Expand Up @@ -346,6 +349,14 @@ EigenProblem::computeResidualAB(const NumericVector<Number> & soln,

computeResidualTags(_fe_vector_tags);

if (bxNormProvided())
{
if (tagA == _nl_eigen->eigenVectorTag())
residualA *= (-*_bx_norm);
if (tagB == _nl_eigen->eigenVectorTag())
residualB *= (-*_bx_norm);
}

_nl_eigen->disassociateVectorFromTag(residualA, tagA);
_nl_eigen->disassociateVectorFromTag(residualB, tagB);
}
Expand All @@ -368,10 +379,7 @@ EigenProblem::computeResidualL2Norm()
_nl_eigen->residualVectorBX() *= eigenvalue;

// Compute entire residual
if (_negative_sign_eigen_kernel)
_nl_eigen->residualVectorAX() += _nl_eigen->residualVectorBX();
else
_nl_eigen->residualVectorAX() -= _nl_eigen->residualVectorBX();
_nl_eigen->residualVectorAX() += _nl_eigen->residualVectorBX();

return _nl_eigen->residualVectorAX().l2_norm();
}
Expand Down Expand Up @@ -642,6 +650,9 @@ EigenProblem::init()
_nl_eigen->sys().use_shell_precond_matrix(solverParams()._precond_matrix_free);
#endif

if (bxNormProvided())
_bx_norm = &getPostprocessorValueByName(*_bx_norm_name);

FEProblemBase::init();
}

Expand Down Expand Up @@ -670,11 +681,11 @@ EigenProblem::initPetscOutputAndSomeSolverSettings()
_app.getOutputWarehouse().solveSetup();
}

Real
EigenProblem::formNorm()
const Real &
EigenProblem::formNorm() const
{
mooseAssert(_bx_norm_name,
"We should not get here unless a bx_norm postprocessor has been provided");
return getPostprocessorValueByName(*_bx_norm_name);
return *_bx_norm;
}
#endif
30 changes: 0 additions & 30 deletions framework/src/utils/SlepcSupport.C
Original file line number Diff line number Diff line change
Expand Up @@ -866,12 +866,6 @@ mooseSlepcEigenFormJacobianB(SNES snes, Vec x, Mat jac, Mat pc, void * ctx)

moosePetscSNESFormMatrixTag(snes, x, pc, sys.get_matrix_B(), ctx, eigen_nl.eigenMatrixTag());

if (eigen_problem->negativeSignEigenKernel())
{
ierr = MatScale(pc, -1.);
CHKERRQ(ierr);
}

PetscFunctionReturn(PETSC_SUCCESS);
}

Expand Down Expand Up @@ -955,12 +949,6 @@ mooseSlepcEigenFormFunctionB(SNES snes, Vec x, Vec r, void * ctx)
else
moosePetscSNESFormFunction(snes, x, r, ctx, eigen_nl.eigenVectorTag());

if (eigen_problem->negativeSignEigenKernel())
{
ierr = VecScale(r, -1.);
CHKERRQ(ierr);
}

PetscFunctionReturn(PETSC_SUCCESS);
}

Expand Down Expand Up @@ -1000,12 +988,6 @@ mooseSlepcEigenFormFunctionAB(SNES /*snes*/, Vec x, Vec Ax, Vec Bx, void * ctx)
ierr = MatMult(B, x, Bx);
CHKERRQ(ierr);

if (eigen_problem->negativeSignEigenKernel())
{
ierr = VecScale(Bx, -1.);
CHKERRQ(ierr);
}

PetscFunctionReturn(PETSC_SUCCESS);
}

Expand All @@ -1030,12 +1012,6 @@ mooseSlepcEigenFormFunctionAB(SNES /*snes*/, Vec x, Vec Ax, Vec Bx, void * ctx)
sys.copy_super_to_sub(*BX, sub_Bx);
}

if (eigen_problem->negativeSignEigenKernel())
{
ierr = VecScale(Bx, -1.);
CHKERRQ(ierr);
}

PetscFunctionReturn(PETSC_SUCCESS);
}

Expand Down Expand Up @@ -1125,12 +1101,6 @@ mooseMatMult_Eigen(Mat mat, Vec x, Vec r)

evaluateResidual(*eigen_problem, x, r, eigen_nl.eigenVectorTag());

if (eigen_problem->negativeSignEigenKernel())
{
ierr = VecScale(r, -1.);
CHKERRQ(ierr);
}

PetscFunctionReturn(PETSC_SUCCESS);
}

Expand Down
49 changes: 41 additions & 8 deletions test/tests/tag/eigen_tag.i
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,10 @@
order = FIRST
family = LAGRANGE
[]
[vec_tag_combined]
order = FIRST
family = LAGRANGE
[]
[mat_tag_diff]
order = FIRST
family = LAGRANGE
Expand All @@ -31,13 +35,14 @@
[diff]
type = Diffusion
variable = u
extra_vector_tags = 'tag_diff'
extra_vector_tags = 'tag_diff tag_combined'
extra_matrix_tags = 'tag_diff'
[]
[rhs]
type = CoefReaction
variable = u
extra_vector_tags = 'eigen tag_rhs'
coefficient = -1
extra_vector_tags = 'eigen tag_rhs tag_combined'
extra_matrix_tags = 'tag_rhs'
[]
[]
Expand All @@ -55,6 +60,12 @@
v = u
vector_tag = tag_rhs
[]
[vec_tag_combined]
type = TagVectorAux
variable = vec_tag_combined
v = u
vector_tag = tag_combined
[]

[mat_tag_diff]
type = TagVectorAux
Expand All @@ -70,18 +81,40 @@
[]
[]

[BCs/homogeneous]
type = DirichletBC
boundary = 'top right bottom left'
variable = u
value = 0
[BCs]
[homogeneous]
type = DirichletBC
boundary = 'top right bottom left'
variable = u
value = 0
extra_vector_tags = 'tag_combined'
[]
[eigen]
type = EigenDirichletBC
variable = u
boundary = 'top right bottom left'
[]
[]

[Problem]
extra_tag_vectors = 'tag_diff tag_rhs'
type = EigenProblem
extra_tag_vectors = 'tag_diff tag_rhs tag_combined'
extra_tag_matrices = 'tag_diff tag_rhs'
bx_norm = uint
[]

[Postprocessors]
[combined_res_norm]
# because all kernels, BCs are tagged, this postprocessor is expected to be zero within tolerance
type = ElementL2Norm
variable = vec_tag_combined
[]
[uint]
type = ElementIntegralVariablePostprocessor
variable = u
execute_on = 'initial linear'
[]
[]

[Executioner]
type = Eigenvalue
Expand Down
Binary file modified test/tests/tag/gold/eigen_tag_out.e
Binary file not shown.
4 changes: 2 additions & 2 deletions test/tests/tag/tests
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
[Tests]
issues = '#9669'
issues = '#9669 #28984'
design = 'interfaces/TaggingInterface.md'

[general]
Expand Down Expand Up @@ -147,7 +147,7 @@
input = 'eigen_tag.i'
exodiff = 'eigen_tag_out.e'
slepc_version = '>=3.8.0'
detail = 'the eigen system, and '
detail = 'the eigen system with and without combined non-eigen and eigen residual objects, and '
[]

[test_old_eigen]
Expand Down

0 comments on commit 7948ad5

Please sign in to comment.