diff --git a/framework/include/kernels/MassEigenKernel.h b/framework/include/kernels/MassEigenKernel.h index 2658598b9e1b..30374cf76080 100644 --- a/framework/include/kernels/MassEigenKernel.h +++ b/framework/include/kernels/MassEigenKernel.h @@ -21,4 +21,5 @@ class MassEigenKernel : public EigenKernel protected: virtual Real computeQpResidual() override; virtual Real computeQpJacobian() override; + const Real _coef; }; diff --git a/framework/include/problems/EigenProblem.h b/framework/include/problems/EigenProblem.h index daca44ba3eb6..bae7c1d7f04e 100644 --- a/framework/include/problems/EigenProblem.h +++ b/framework/include/problems/EigenProblem.h @@ -224,21 +224,15 @@ class EigenProblem : public FEProblemBase */ void wereMatricesFormed(bool mf) { _matrices_formed = mf; } -private: /** - * Do some free/extra power iterations + * Form the Bx norm */ - void doFreeNonlinearPowerIterations(unsigned int free_power_iterations); + Real formNorm(); /** - * Adjust eigen vector by either scaling the existing values or setting new values - * The operations are applied for only eigen variables + * Whether a Bx norm postprocessor has been provided */ - void adjustEigenVector(const Real value, bool scaling); - -#endif - - using FEProblemBase::_nl; + bool bxNormProvided() const { return _bx_norm_name.has_value(); } protected: unsigned int _n_eigen_pairs_required; @@ -275,6 +269,26 @@ class EigenProblem : public FEProblemBase bool & _first_solve; /// A value used for initial normalization Real _initial_eigenvalue; + +private: + /** + * Do some free/extra power iterations + */ + void doFreeNonlinearPowerIterations(unsigned int free_power_iterations); + + /** + * Adjust eigen vector by either scaling the existing values or setting new values + * The operations are applied for only eigen variables + */ + void adjustEigenVector(const Real value, bool scaling); + + /// The name of the Postprocessor providing the Bx norm. This may be empty in which case the + /// default L2 norm of Bx will be used as the Bx norm + const std::optional _bx_norm_name; + +#endif + + using FEProblemBase::_nl; }; #ifdef LIBMESH_HAVE_SLEPC diff --git a/framework/src/kernels/MassEigenKernel.C b/framework/src/kernels/MassEigenKernel.C index 1ec02515e8ba..dd5006040eaa 100644 --- a/framework/src/kernels/MassEigenKernel.C +++ b/framework/src/kernels/MassEigenKernel.C @@ -15,21 +15,26 @@ InputParameters MassEigenKernel::validParams() { InputParameters params = EigenKernel::validParams(); - params.addClassDescription("An eigenkernel with weak form $\\lambda(\\psi_i, -u_h)$ where " - "$\\lambda$ is the eigenvalue."); + params.addClassDescription( + "An eigenkernel with weak form $\\lambda(\\psi_i, -u_h \\times coeff)$ where " + "$\\lambda$ is the eigenvalue."); + params.addParam("coefficient", 1.0, "Coefficient multiplying the term"); return params; } -MassEigenKernel::MassEigenKernel(const InputParameters & parameters) : EigenKernel(parameters) {} +MassEigenKernel::MassEigenKernel(const InputParameters & parameters) + : EigenKernel(parameters), _coef(this->template getParam("coefficient")) +{ +} Real MassEigenKernel::computeQpResidual() { - return -_u[_qp] * _test[_i][_qp]; + return -_coef * _u[_qp] * _test[_i][_qp]; } Real MassEigenKernel::computeQpJacobian() { - return -_phi[_j][_qp] * _test[_i][_qp]; + return -_coef * _phi[_j][_qp] * _test[_i][_qp]; } diff --git a/framework/src/problems/EigenProblem.C b/framework/src/problems/EigenProblem.C index 5d0678eb3d17..cfce0f6d9d81 100644 --- a/framework/src/problems/EigenProblem.C +++ b/framework/src/problems/EigenProblem.C @@ -42,6 +42,7 @@ EigenProblem::validParams() "active_eigen_index", 0, "Which eigenvector is used to compute residual and also associated to nonlinear variable"); + params.addParam("bx_norm", "A postprocessor describing the norm of Bx"); return params; } @@ -60,7 +61,10 @@ EigenProblem::EigenProblem(const InputParameters & parameters) _constant_matrices(false), _has_normalization(false), _normal_factor(1.0), - _first_solve(declareRestartableData("first_solve", true)) + _first_solve(declareRestartableData("first_solve", true)), + _bx_norm_name(isParamValid("bx_norm") + ? std::make_optional(getParam("bx_norm")) + : std::nullopt) { #ifdef LIBMESH_HAVE_SLEPC if (_nl_sys_names.size() > 1) @@ -416,7 +420,7 @@ EigenProblem::preScaleEigenVector(const std::pair & eig) // Eigenvalue magnitude Real v = std::sqrt(eig.first * eig.first + eig.second * eig.second); // Scaling factor - Real factor = 1 / v / _nl_eigen->residualVectorBX().l2_norm(); + Real factor = 1 / v / (bxNormProvided() ? formNorm() : _nl_eigen->residualVectorBX().l2_norm()); // Scale eigenvector if (!MooseUtils::absoluteFuzzyEqual(factor, 1)) scaleEigenvector(factor); @@ -637,4 +641,11 @@ EigenProblem::initPetscOutput() _app.getOutputWarehouse().solveSetup(); } +Real +EigenProblem::formNorm() +{ + mooseAssert(_bx_norm_name, + "We should not get here unless a bx_norm postprocessor has been provided"); + return getPostprocessorValueByName(*_bx_norm_name); +} #endif diff --git a/framework/src/utils/SlepcSupport.C b/framework/src/utils/SlepcSupport.C index 10fa6fc48445..a9d54b0e40ac 100644 --- a/framework/src/utils/SlepcSupport.C +++ b/framework/src/utils/SlepcSupport.C @@ -356,6 +356,9 @@ setFreeNonlinearPowerIterations(unsigned int free_power_iterations) // -snes_no_convergence_test is a perfect option, but it was removed from PETSc Moose::PetscSupport::setSinglePetscOption("-snes_rtol", "0.99999999999"); Moose::PetscSupport::setSinglePetscOption("-eps_max_it", stringify(free_power_iterations)); + // We always want the number of free power iterations respected so we don't want to stop early if + // we've satisfied a convergence criterion. Consequently we make this tolerance very tight + Moose::PetscSupport::setSinglePetscOption("-eps_tol", "1e-50"); } void @@ -367,6 +370,7 @@ clearFreeNonlinearPowerIterations(const InputParameters & params) stringify(params.get("nl_max_its"))); Moose::PetscSupport::setSinglePetscOption("-snes_rtol", stringify(params.get("nl_rel_tol"))); + Moose::PetscSupport::setSinglePetscOption("-eps_tol", stringify(params.get("eigen_tol"))); } void @@ -962,28 +966,60 @@ mooseSlepcEigenFormFunctionAB(SNES /*snes*/, Vec x, Vec Ax, Vec Bx, void * ctx) PetscFunctionReturn(0); } +PetscErrorCode +mooseSlepcEigenFormNorm(SNES /*snes*/, Vec /*Bx*/, PetscReal * norm, void * ctx) +{ + PetscFunctionBegin; + auto * const eigen_problem = static_cast(ctx); + *norm = eigen_problem->formNorm(); + PetscFunctionReturn(0); +} + void attachCallbacksToMat(EigenProblem & eigen_problem, Mat mat, bool eigen) { + // Recall that we are solving the potentially nonlinear problem: + // F(x) = A(x) - \lambda B(x) = 0 + // + // To solve this, we can use Newton's method: J \Delta x = -F + // Generally we will approximate J using matrix free methods. However, in order to solve the + // linearized system efficiently, we typically will need preconditioning. Typically we will build + // the preconditioner only from A, but we also have the option to include information from B + + // Attach the Jacobian computation function. If \p mat is the "eigen" matrix corresponding to B, + // then attach our JacobianB computation routine, else the matrix corresponds to A, and we attach + // the JacobianA computation routine PetscObjectComposeFunction((PetscObject)mat, "formJacobian", eigen ? Moose::SlepcSupport::mooseSlepcEigenFormJacobianB : Moose::SlepcSupport::mooseSlepcEigenFormJacobianA); + + // Attach the residual computation function. If \p mat is the "eigen" matrix corresponding to B, + // then attach our FunctionB computation routine, else the matrix corresponds to A, and we attach + // the FunctionA computation routine PetscObjectComposeFunction((PetscObject)mat, "formFunction", eigen ? Moose::SlepcSupport::mooseSlepcEigenFormFunctionB : Moose::SlepcSupport::mooseSlepcEigenFormFunctionA); + // It's also beneficial to be able to evaluate both A and B residuals at once PetscObjectComposeFunction( (PetscObject)mat, "formFunctionAB", Moose::SlepcSupport::mooseSlepcEigenFormFunctionAB); + // Users may choose to provide a custom measure of the norm of B (Bx for a linear system) + if (eigen_problem.bxNormProvided()) + PetscObjectComposeFunction( + (PetscObject)mat, "formNorm", Moose::SlepcSupport::mooseSlepcEigenFormNorm); + + // Finally we need to attach the "context" object, which is our EigenProblem, to the matrices so + // that eventually when we get callbacks from SLEPc we can call methods on the EigenProblem PetscContainer container; PetscContainerCreate(eigen_problem.comm().get(), &container); PetscContainerSetPointer(container, &eigen_problem); - PetscObjectCompose((PetscObject)mat, "formJacobianCtx", nullptr); PetscObjectCompose((PetscObject)mat, "formJacobianCtx", (PetscObject)container); - PetscObjectCompose((PetscObject)mat, "formFunctionCtx", nullptr); PetscObjectCompose((PetscObject)mat, "formFunctionCtx", (PetscObject)container); + if (eigen_problem.bxNormProvided()) + PetscObjectCompose((PetscObject)mat, "formNormCtx", (PetscObject)container); PetscContainerDestroy(&container); } @@ -1086,7 +1122,8 @@ PCCreate_MoosePC(PC pc) PetscFunctionReturn(0); } -PetscErrorCode PCDestroy_MoosePC(PC /*pc*/) +PetscErrorCode +PCDestroy_MoosePC(PC /*pc*/) { PetscFunctionBegin; /* We do not need to do anything right now, but later we may have some data we need to free here diff --git a/test/tests/executioners/eigen_convergence/a.i b/test/tests/executioners/eigen_convergence/a.i new file mode 100644 index 000000000000..69a1aad33486 --- /dev/null +++ b/test/tests/executioners/eigen_convergence/a.i @@ -0,0 +1,84 @@ +[Mesh] + [gen] + type = GeneratedMeshGenerator + dim = 2 + xmin = 0 + xmax = 160 + ymin = 0 + ymax = 160 + nx = 8 + ny = 8 + [] + uniform_refine = 0 +[] + +[Variables] + [u] + [] +[] + +[Kernels] + [diffusion] + type = MatDiffusion + diffusivity = diffusivity + variable = u + [] + [reaction] + type = CoefReaction + coefficient = 0.01 + variable = u + [] + [rhs] + type = CoefReaction + extra_vector_tags = 'eigen' + coefficient = -0.01 + variable = u + [] +[] + +[BCs] + [robin] + type = VacuumBC + boundary = 'left bottom' + variable = u + [] +[] + +[Materials] + [nm] + type = GenericConstantMaterial + block = 0 + prop_names = 'diffusivity' + prop_values = 0.333333333333333333 + [] +[] + +[VectorPostprocessors] + [eigen] + type = Eigenvalues + inverse_eigenvalue = true + [] +[] + +[Postprocessors] + [fluxintegral] + type = ElementIntegralVariablePostprocessor + variable = u + execute_on = linear + [] +[] + +[Problem] + type = EigenProblem +[] + +[Executioner] + type = Eigenvalue + solve_type = PJFNK + free_power_iterations = 4 + nl_abs_tol = 2e-10 +[] + +[Outputs] + csv = true +[] diff --git a/test/tests/executioners/eigen_convergence/b.i b/test/tests/executioners/eigen_convergence/b.i new file mode 100644 index 000000000000..26b4bbfa87f3 --- /dev/null +++ b/test/tests/executioners/eigen_convergence/b.i @@ -0,0 +1,73 @@ +[Mesh] + [gen] + type = GeneratedMeshGenerator + dim = 2 + xmin = 0 + xmax = 160 + ymin = 0 + ymax = 160 + nx = 8 + ny = 8 + [] + uniform_refine = 0 +[] + +[Variables] + [u] + [] +[] + +[Kernels] + [diffusion] + type = MatDiffusion + diffusivity = diffusivity + variable = u + [] + [reaction] + type = CoefReaction + coefficient = 0.01 + variable = u + [] + [rhs] + type = MassEigenKernel + coefficient = 0.01 + variable = u + [] +[] + +[BCs] + [robin] + type = VacuumBC + boundary = 'left bottom' + variable = u + [] +[] + +[Materials] + [nm] + type = GenericConstantMaterial + block = 0 + prop_names = 'diffusivity' + prop_values = 0.333333333333333333 + [] +[] + +[Postprocessors] + [fluxintegral] + type = ElementIntegralVariablePostprocessor + variable = u + execute_on = linear + [] +[] + +[Executioner] + type = NonlinearEigen + bx_norm = fluxintegral + solve_type = PJFNK + free_power_iterations = 4 + nl_abs_tol = 2e-10 +[] + +[Outputs] + csv = true +[] diff --git a/test/tests/executioners/eigen_convergence/gold/a_out_eigen_0001.csv b/test/tests/executioners/eigen_convergence/gold/a_out_eigen_0001.csv new file mode 100644 index 000000000000..df5723e7c3b5 --- /dev/null +++ b/test/tests/executioners/eigen_convergence/gold/a_out_eigen_0001.csv @@ -0,0 +1,2 @@ +eigen_values_imag,eigen_values_real +0,0.99364796564363 diff --git a/test/tests/executioners/eigen_convergence/gold/b_out.csv b/test/tests/executioners/eigen_convergence/gold/b_out.csv new file mode 100644 index 000000000000..a2649e3df945 --- /dev/null +++ b/test/tests/executioners/eigen_convergence/gold/b_out.csv @@ -0,0 +1,4 @@ +time,fluxintegral,Executioner/eigenvalue +0,1,1 +1,0.97869657767363,0.97869657767363 +2,0.99364796612084,0.99364796612084 diff --git a/test/tests/executioners/eigen_convergence/tests b/test/tests/executioners/eigen_convergence/tests new file mode 100644 index 000000000000..1848c1dd43e0 --- /dev/null +++ b/test/tests/executioners/eigen_convergence/tests @@ -0,0 +1,43 @@ +[Tests] + design = 'EigenProblem.md' + issues = '#20095' + [testb] + type = CSVDiff + csvdiff = 'b_out.csv' + input = 'b.i' + expect_out = '4 Nonlinear' + absent_out = '5 Nonlinear' + requirement = 'The system shall be able to solve an eigenvalue problem using something other than the L2 norm of Bx for normalization and a native eigenvalue solver.' + rel_err = 3e-5 + [] + [test11] + type = CSVDiff + csvdiff = 'a_out_eigen_0001.csv' + input = 'a.i' + expect_out = '1[01] Nonlinear' + absent_out = '12 Nonlinear' + requirement = 'The system shall be able to solve an eigenvalue problem using the L2 norm of Bx and the sign of its first nonzero entry for normalization with a SLEPc eigenvalue solver.' + [] + [test7] + type = CSVDiff + csvdiff = 'a_out_eigen_0001.csv' + prereq = test11 + input = 'a.i' + expect_out = '7 Nonlinear' + absent_out = '8 Nonlinear' + cli_args = 'Problem/bx_norm=fluxintegral' + petsc_version = '>=3.20.0' + requirement = 'The system shall be able to solve an eigenvalue problem using the sign of the first nonzero entry of Bx combined with something other than the L2 norm of Bx for normalization with a SLEPc eigenvalue solver.' + [] + [test4] + type = CSVDiff + csvdiff = 'a_out_eigen_0001.csv' + prereq = test7 + input = 'a.i' + expect_out = '4 Nonlinear' + absent_out = '5 Nonlinear' + cli_args = 'Problem/bx_norm=fluxintegral -eps_power_sign_normalization 0' + petsc_version = '>=3.20.0' + requirement = 'The system shall be able to solve an eigenvalue problem using something other than the L2 norm of Bx for normalization and a SLEPc eigenvalue solver.' + [] +[]