From 2cfe5b3223ffeabbf3981fe8c48e027df22e886b Mon Sep 17 00:00:00 2001 From: Yaqi Wang Date: Wed, 12 Oct 2022 13:27:43 -0500 Subject: [PATCH 01/16] demonstration of Issue #20095 --- framework/include/kernels/MassEigenKernel.h | 1 + framework/src/kernels/MassEigenKernel.C | 10 ++- .../executioners/eigen_convergence_issue/a.i | 74 +++++++++++++++++++ .../executioners/eigen_convergence_issue/b.i | 74 +++++++++++++++++++ 4 files changed, 156 insertions(+), 3 deletions(-) create mode 100644 test/tests/executioners/eigen_convergence_issue/a.i create mode 100644 test/tests/executioners/eigen_convergence_issue/b.i 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/src/kernels/MassEigenKernel.C b/framework/src/kernels/MassEigenKernel.C index 1ec02515e8ba..0151a31c9e2c 100644 --- a/framework/src/kernels/MassEigenKernel.C +++ b/framework/src/kernels/MassEigenKernel.C @@ -17,19 +17,23 @@ MassEigenKernel::validParams() InputParameters params = EigenKernel::validParams(); params.addClassDescription("An eigenkernel with weak form $\\lambda(\\psi_i, -u_h)$ where " "$\\lambda$ is the eigenvalue."); + params.addParam("coefficient", 1.0, "Coefficient of 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/test/tests/executioners/eigen_convergence_issue/a.i b/test/tests/executioners/eigen_convergence_issue/a.i new file mode 100644 index 000000000000..de5e10ce24fb --- /dev/null +++ b/test/tests/executioners/eigen_convergence_issue/a.i @@ -0,0 +1,74 @@ +[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 + [] +[] + +[Postprocessors] + [fluxintegral] + type = ElementIntegralVariablePostprocessor + variable = u + execute_on = linear + [] +[] + +[Executioner] + type = Eigenvalue + solve_type = PJFNK + free_power_iterations = 4 + nl_abs_tol = 1e-10 +[] + +[Outputs] + exodus = true + perf_graph = true +[] diff --git a/test/tests/executioners/eigen_convergence_issue/b.i b/test/tests/executioners/eigen_convergence_issue/b.i new file mode 100644 index 000000000000..b1c4332e638b --- /dev/null +++ b/test/tests/executioners/eigen_convergence_issue/b.i @@ -0,0 +1,74 @@ +[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 = 1e-10 +[] + +[Outputs] + exodus = true + perf_graph = true +[] From 7bf542a459e62ad2b5e0c4c89840945a9b544101 Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Mon, 24 Jul 2023 17:44:28 -0700 Subject: [PATCH 02/16] Testing new formNorm from SLEPc. Getting incorrect eigenvalue --- framework/include/problems/EigenProblem.h | 33 +++++++++++++------ framework/src/problems/EigenProblem.C | 11 ++++++- framework/src/utils/SlepcSupport.C | 17 ++++++++++ .../executioners/eigen_convergence_issue/a.i | 5 +++ 4 files changed, 55 insertions(+), 11 deletions(-) diff --git a/framework/include/problems/EigenProblem.h b/framework/include/problems/EigenProblem.h index daca44ba3eb6..6fe43dd8f456 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_provided; } protected: unsigned int _n_eigen_pairs_required; @@ -275,6 +269,25 @@ 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); + + /// Whether the bx norm postprocessor is provided + const bool _bx_norm_provided; + +#endif + + using FEProblemBase::_nl; }; #ifdef LIBMESH_HAVE_SLEPC diff --git a/framework/src/problems/EigenProblem.C b/framework/src/problems/EigenProblem.C index 5d0678eb3d17..6af5b86500eb 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,8 @@ 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_provided(isParamValid("bx_norm")) { #ifdef LIBMESH_HAVE_SLEPC if (_nl_sys_names.size() > 1) @@ -637,4 +639,11 @@ EigenProblem::initPetscOutput() _app.getOutputWarehouse().solveSetup(); } +Real +EigenProblem::formNorm() +{ + mooseAssert(_bx_norm_provided, + "We should not get here unless a bx_norm postprocessor has been provided"); + return getPostprocessorValueByName(getParam("bx_norm")); +} #endif diff --git a/framework/src/utils/SlepcSupport.C b/framework/src/utils/SlepcSupport.C index 10fa6fc48445..6c2fe0b8ad5f 100644 --- a/framework/src/utils/SlepcSupport.C +++ b/framework/src/utils/SlepcSupport.C @@ -962,6 +962,15 @@ 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(PETSC_SUCCESS); +} + void attachCallbacksToMat(EigenProblem & eigen_problem, Mat mat, bool eigen) { @@ -976,6 +985,9 @@ attachCallbacksToMat(EigenProblem & eigen_problem, Mat mat, bool eigen) PetscObjectComposeFunction( (PetscObject)mat, "formFunctionAB", Moose::SlepcSupport::mooseSlepcEigenFormFunctionAB); + if (eigen_problem.bxNormProvided()) + PetscObjectComposeFunction( + (PetscObject)mat, "formNorm", Moose::SlepcSupport::mooseSlepcEigenFormNorm); PetscContainer container; PetscContainerCreate(eigen_problem.comm().get(), &container); @@ -984,6 +996,11 @@ attachCallbacksToMat(EigenProblem & eigen_problem, Mat mat, bool eigen) 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", nullptr); + PetscObjectCompose((PetscObject)mat, "formNormCtx", (PetscObject)container); + } PetscContainerDestroy(&container); } diff --git a/test/tests/executioners/eigen_convergence_issue/a.i b/test/tests/executioners/eigen_convergence_issue/a.i index de5e10ce24fb..b55a5053549a 100644 --- a/test/tests/executioners/eigen_convergence_issue/a.i +++ b/test/tests/executioners/eigen_convergence_issue/a.i @@ -61,6 +61,11 @@ [] [] +[Problem] + type = EigenProblem + bx_norm = fluxintegral +[] + [Executioner] type = Eigenvalue solve_type = PJFNK From 67a00aa93b2caf92af36e9c04b0203db1004ab8b Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Tue, 25 Jul 2023 17:18:22 -0700 Subject: [PATCH 03/16] Use our custom bx_norm for pre-scaling --- framework/src/problems/EigenProblem.C | 2 +- test/tests/executioners/eigen_convergence_issue/a.i | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/framework/src/problems/EigenProblem.C b/framework/src/problems/EigenProblem.C index 6af5b86500eb..affde4dd2ad8 100644 --- a/framework/src/problems/EigenProblem.C +++ b/framework/src/problems/EigenProblem.C @@ -418,7 +418,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 / (_bx_norm_provided ? formNorm() : _nl_eigen->residualVectorBX().l2_norm()); // Scale eigenvector if (!MooseUtils::absoluteFuzzyEqual(factor, 1)) scaleEigenvector(factor); diff --git a/test/tests/executioners/eigen_convergence_issue/a.i b/test/tests/executioners/eigen_convergence_issue/a.i index b55a5053549a..bf6a5c2283a0 100644 --- a/test/tests/executioners/eigen_convergence_issue/a.i +++ b/test/tests/executioners/eigen_convergence_issue/a.i @@ -71,6 +71,8 @@ solve_type = PJFNK free_power_iterations = 4 nl_abs_tol = 1e-10 + petsc_options_iname = '-eps_tol' + petsc_options_value = '1e-10' [] [Outputs] From fef33cf89947fb2c0272409158cb5bb8bd461220 Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Thu, 3 Aug 2023 08:46:11 -0700 Subject: [PATCH 04/16] Add tests for eigen nonlinear convergence --- .../executioners/eigen_convergence_issue/a.i | 6 ++-- .../eigen_convergence_issue/tests | 29 +++++++++++++++++++ 2 files changed, 32 insertions(+), 3 deletions(-) create mode 100644 test/tests/executioners/eigen_convergence_issue/tests diff --git a/test/tests/executioners/eigen_convergence_issue/a.i b/test/tests/executioners/eigen_convergence_issue/a.i index bf6a5c2283a0..29808363e8a7 100644 --- a/test/tests/executioners/eigen_convergence_issue/a.i +++ b/test/tests/executioners/eigen_convergence_issue/a.i @@ -63,7 +63,6 @@ [Problem] type = EigenProblem - bx_norm = fluxintegral [] [Executioner] @@ -71,8 +70,9 @@ solve_type = PJFNK free_power_iterations = 4 nl_abs_tol = 1e-10 - petsc_options_iname = '-eps_tol' - petsc_options_value = '1e-10' + # The first set of power iterations don't seem to respect these options + # petsc_options_iname = '-eps_tol' + # petsc_options_value = '1e-10' [] [Outputs] diff --git a/test/tests/executioners/eigen_convergence_issue/tests b/test/tests/executioners/eigen_convergence_issue/tests new file mode 100644 index 000000000000..b7536f364c58 --- /dev/null +++ b/test/tests/executioners/eigen_convergence_issue/tests @@ -0,0 +1,29 @@ +[Tests] + [testb] + type = RunApp + input = 'b.i' + expect_out = '4 Nonlinear' + absent_out = '5 Nonlinear' + [] + [test11] + type = RunApp + input = 'a.i' + expect_out = '11 Nonlinear' + absent_out = '12 Nonlinear' + cli_args = '-eps_tol 1e-10' + [] + [test7] + type = RunApp + input = 'a.i' + expect_out = '7 Nonlinear' + absent_out = '8 Nonlinear' + cli_args = 'Problem/bx_norm=fluxintegral -eps_tol 1e-10' + [] + [test4] + type = RunApp + input = 'a.i' + expect_out = '4 Nonlinear' + absent_out = '5 Nonlinear' + cli_args = 'Problem/bx_norm=fluxintegral -eps_power_sign_normalization 0 -eps_tol 1e-10' + [] +[] From 1b0fd0c21b39b99da93fd6aa3a45fc67d0f3aeba Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Thu, 2 Nov 2023 14:13:41 -0700 Subject: [PATCH 05/16] Adjust tolerance for robustness --- test/tests/executioners/eigen_convergence_issue/a.i | 2 +- test/tests/executioners/eigen_convergence_issue/b.i | 2 +- test/tests/executioners/eigen_convergence_issue/tests | 6 +++--- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/test/tests/executioners/eigen_convergence_issue/a.i b/test/tests/executioners/eigen_convergence_issue/a.i index 29808363e8a7..80d41e95bc03 100644 --- a/test/tests/executioners/eigen_convergence_issue/a.i +++ b/test/tests/executioners/eigen_convergence_issue/a.i @@ -69,7 +69,7 @@ type = Eigenvalue solve_type = PJFNK free_power_iterations = 4 - nl_abs_tol = 1e-10 + nl_abs_tol = 2e-10 # The first set of power iterations don't seem to respect these options # petsc_options_iname = '-eps_tol' # petsc_options_value = '1e-10' diff --git a/test/tests/executioners/eigen_convergence_issue/b.i b/test/tests/executioners/eigen_convergence_issue/b.i index b1c4332e638b..edaa4a87534e 100644 --- a/test/tests/executioners/eigen_convergence_issue/b.i +++ b/test/tests/executioners/eigen_convergence_issue/b.i @@ -65,7 +65,7 @@ bx_norm = fluxintegral solve_type = PJFNK free_power_iterations = 4 - nl_abs_tol = 1e-10 + nl_abs_tol = 2e-10 [] [Outputs] diff --git a/test/tests/executioners/eigen_convergence_issue/tests b/test/tests/executioners/eigen_convergence_issue/tests index b7536f364c58..77f591a979c7 100644 --- a/test/tests/executioners/eigen_convergence_issue/tests +++ b/test/tests/executioners/eigen_convergence_issue/tests @@ -10,20 +10,20 @@ input = 'a.i' expect_out = '11 Nonlinear' absent_out = '12 Nonlinear' - cli_args = '-eps_tol 1e-10' + cli_args = '-eps_tol 2e-10' [] [test7] type = RunApp input = 'a.i' expect_out = '7 Nonlinear' absent_out = '8 Nonlinear' - cli_args = 'Problem/bx_norm=fluxintegral -eps_tol 1e-10' + cli_args = 'Problem/bx_norm=fluxintegral -eps_tol 2e-10' [] [test4] type = RunApp input = 'a.i' expect_out = '4 Nonlinear' absent_out = '5 Nonlinear' - cli_args = 'Problem/bx_norm=fluxintegral -eps_power_sign_normalization 0 -eps_tol 1e-10' + cli_args = 'Problem/bx_norm=fluxintegral -eps_power_sign_normalization 0 -eps_tol 2e-10' [] [] From 5ebc728546ad47f4248e5fc527a62422855faf2f Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Thu, 2 Nov 2023 15:48:10 -0700 Subject: [PATCH 06/16] Protect against old petsc versions --- framework/src/utils/SlepcSupport.C | 5 +++-- test/tests/executioners/eigen_convergence_issue/tests | 2 ++ 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/framework/src/utils/SlepcSupport.C b/framework/src/utils/SlepcSupport.C index 6c2fe0b8ad5f..9e2122e61be9 100644 --- a/framework/src/utils/SlepcSupport.C +++ b/framework/src/utils/SlepcSupport.C @@ -968,7 +968,7 @@ mooseSlepcEigenFormNorm(SNES /*snes*/, Vec /*Bx*/, PetscReal * norm, void * ctx) PetscFunctionBegin; auto * const eigen_problem = static_cast(ctx); *norm = eigen_problem->formNorm(); - PetscFunctionReturn(PETSC_SUCCESS); + PetscFunctionReturn(0); } void @@ -1103,7 +1103,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_issue/tests b/test/tests/executioners/eigen_convergence_issue/tests index 77f591a979c7..0d63b4a516ee 100644 --- a/test/tests/executioners/eigen_convergence_issue/tests +++ b/test/tests/executioners/eigen_convergence_issue/tests @@ -18,6 +18,7 @@ expect_out = '7 Nonlinear' absent_out = '8 Nonlinear' cli_args = 'Problem/bx_norm=fluxintegral -eps_tol 2e-10' + petsc_version = '>=3.20.0' [] [test4] type = RunApp @@ -25,5 +26,6 @@ expect_out = '4 Nonlinear' absent_out = '5 Nonlinear' cli_args = 'Problem/bx_norm=fluxintegral -eps_power_sign_normalization 0 -eps_tol 2e-10' + petsc_version = '>=3.20.0' [] [] From 7949594fef18ba30034ff15e05c74631d1806778 Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Fri, 3 Nov 2023 17:18:45 -0700 Subject: [PATCH 07/16] Some more test edits - Use CSVDiff so we can be confident the eigenvalue is good - Allow 10 or 11 nonlinear iterations for 'test11' - Add design, issues, and requirements --- .../executioners/eigen_convergence_issue/a.i | 10 +++++++-- .../executioners/eigen_convergence_issue/b.i | 3 +-- .../gold/a_out_eigen_0001.csv | 2 ++ .../eigen_convergence_issue/gold/b_out.csv | 4 ++++ .../eigen_convergence_issue/tests | 22 ++++++++++++++----- 5 files changed, 32 insertions(+), 9 deletions(-) create mode 100644 test/tests/executioners/eigen_convergence_issue/gold/a_out_eigen_0001.csv create mode 100644 test/tests/executioners/eigen_convergence_issue/gold/b_out.csv diff --git a/test/tests/executioners/eigen_convergence_issue/a.i b/test/tests/executioners/eigen_convergence_issue/a.i index 80d41e95bc03..f67516a5feb3 100644 --- a/test/tests/executioners/eigen_convergence_issue/a.i +++ b/test/tests/executioners/eigen_convergence_issue/a.i @@ -53,6 +53,13 @@ [] [] +[VectorPostprocessors] + [eigen] + type = Eigenvalues + inverse_eigenvalue = true + [] +[] + [Postprocessors] [fluxintegral] type = ElementIntegralVariablePostprocessor @@ -76,6 +83,5 @@ [] [Outputs] - exodus = true - perf_graph = true + csv = true [] diff --git a/test/tests/executioners/eigen_convergence_issue/b.i b/test/tests/executioners/eigen_convergence_issue/b.i index edaa4a87534e..26b4bbfa87f3 100644 --- a/test/tests/executioners/eigen_convergence_issue/b.i +++ b/test/tests/executioners/eigen_convergence_issue/b.i @@ -69,6 +69,5 @@ [] [Outputs] - exodus = true - perf_graph = true + csv = true [] diff --git a/test/tests/executioners/eigen_convergence_issue/gold/a_out_eigen_0001.csv b/test/tests/executioners/eigen_convergence_issue/gold/a_out_eigen_0001.csv new file mode 100644 index 000000000000..df5723e7c3b5 --- /dev/null +++ b/test/tests/executioners/eigen_convergence_issue/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_issue/gold/b_out.csv b/test/tests/executioners/eigen_convergence_issue/gold/b_out.csv new file mode 100644 index 000000000000..a2649e3df945 --- /dev/null +++ b/test/tests/executioners/eigen_convergence_issue/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_issue/tests b/test/tests/executioners/eigen_convergence_issue/tests index 0d63b4a516ee..dd8fd5d337cf 100644 --- a/test/tests/executioners/eigen_convergence_issue/tests +++ b/test/tests/executioners/eigen_convergence_issue/tests @@ -1,31 +1,43 @@ [Tests] + design = 'EigenProblem.md' + issues = '#20095' [testb] - type = RunApp + 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.' [] [test11] - type = RunApp + type = CSVDiff + csvdiff = 'a_out_eigen_0001.csv' input = 'a.i' - expect_out = '11 Nonlinear' + expect_out = '1[01] Nonlinear' absent_out = '12 Nonlinear' cli_args = '-eps_tol 2e-10' + 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 = RunApp + 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 -eps_tol 2e-10' 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 = RunApp + 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 -eps_tol 2e-10' 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.' [] [] From 8d8b8fcac4d845e2c1be5ece510b0e35be0448de Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Fri, 3 Nov 2023 17:20:25 -0700 Subject: [PATCH 08/16] git mv eigen_convergence_issue eigen_convergence --- .../{eigen_convergence_issue => eigen_convergence}/a.i | 0 .../{eigen_convergence_issue => eigen_convergence}/b.i | 0 .../gold/a_out_eigen_0001.csv | 0 .../{eigen_convergence_issue => eigen_convergence}/gold/b_out.csv | 0 .../{eigen_convergence_issue => eigen_convergence}/tests | 0 5 files changed, 0 insertions(+), 0 deletions(-) rename test/tests/executioners/{eigen_convergence_issue => eigen_convergence}/a.i (100%) rename test/tests/executioners/{eigen_convergence_issue => eigen_convergence}/b.i (100%) rename test/tests/executioners/{eigen_convergence_issue => eigen_convergence}/gold/a_out_eigen_0001.csv (100%) rename test/tests/executioners/{eigen_convergence_issue => eigen_convergence}/gold/b_out.csv (100%) rename test/tests/executioners/{eigen_convergence_issue => eigen_convergence}/tests (100%) diff --git a/test/tests/executioners/eigen_convergence_issue/a.i b/test/tests/executioners/eigen_convergence/a.i similarity index 100% rename from test/tests/executioners/eigen_convergence_issue/a.i rename to test/tests/executioners/eigen_convergence/a.i diff --git a/test/tests/executioners/eigen_convergence_issue/b.i b/test/tests/executioners/eigen_convergence/b.i similarity index 100% rename from test/tests/executioners/eigen_convergence_issue/b.i rename to test/tests/executioners/eigen_convergence/b.i diff --git a/test/tests/executioners/eigen_convergence_issue/gold/a_out_eigen_0001.csv b/test/tests/executioners/eigen_convergence/gold/a_out_eigen_0001.csv similarity index 100% rename from test/tests/executioners/eigen_convergence_issue/gold/a_out_eigen_0001.csv rename to test/tests/executioners/eigen_convergence/gold/a_out_eigen_0001.csv diff --git a/test/tests/executioners/eigen_convergence_issue/gold/b_out.csv b/test/tests/executioners/eigen_convergence/gold/b_out.csv similarity index 100% rename from test/tests/executioners/eigen_convergence_issue/gold/b_out.csv rename to test/tests/executioners/eigen_convergence/gold/b_out.csv diff --git a/test/tests/executioners/eigen_convergence_issue/tests b/test/tests/executioners/eigen_convergence/tests similarity index 100% rename from test/tests/executioners/eigen_convergence_issue/tests rename to test/tests/executioners/eigen_convergence/tests From 25d08d1e71699e095436280592d842969e690d82 Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Fri, 3 Nov 2023 18:58:14 -0700 Subject: [PATCH 09/16] Loosen tol for csvdiff for legacy eigen --- test/tests/executioners/eigen_convergence/tests | 1 + 1 file changed, 1 insertion(+) diff --git a/test/tests/executioners/eigen_convergence/tests b/test/tests/executioners/eigen_convergence/tests index dd8fd5d337cf..8d0b78d062ef 100644 --- a/test/tests/executioners/eigen_convergence/tests +++ b/test/tests/executioners/eigen_convergence/tests @@ -8,6 +8,7 @@ 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 From ff35e4f9d2e8963502faf33e311d8ad2814a6ef4 Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Wed, 8 Nov 2023 14:52:44 -0700 Subject: [PATCH 10/16] Get name of bx norm pp in constructor --- framework/include/problems/EigenProblem.h | 7 ++++--- framework/src/problems/EigenProblem.C | 8 ++++---- 2 files changed, 8 insertions(+), 7 deletions(-) diff --git a/framework/include/problems/EigenProblem.h b/framework/include/problems/EigenProblem.h index 6fe43dd8f456..d0774674cad4 100644 --- a/framework/include/problems/EigenProblem.h +++ b/framework/include/problems/EigenProblem.h @@ -232,7 +232,7 @@ class EigenProblem : public FEProblemBase /** * Whether a Bx norm postprocessor has been provided */ - bool bxNormProvided() const { return _bx_norm_provided; } + bool bxNormProvided() const { return !_bx_norm_name.empty(); } protected: unsigned int _n_eigen_pairs_required; @@ -282,8 +282,9 @@ class EigenProblem : public FEProblemBase */ void adjustEigenVector(const Real value, bool scaling); - /// Whether the bx norm postprocessor is provided - const bool _bx_norm_provided; + /// 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 PostprocessorName _bx_norm_name; #endif diff --git a/framework/src/problems/EigenProblem.C b/framework/src/problems/EigenProblem.C index affde4dd2ad8..3348bf97a40c 100644 --- a/framework/src/problems/EigenProblem.C +++ b/framework/src/problems/EigenProblem.C @@ -62,7 +62,7 @@ EigenProblem::EigenProblem(const InputParameters & parameters) _has_normalization(false), _normal_factor(1.0), _first_solve(declareRestartableData("first_solve", true)), - _bx_norm_provided(isParamValid("bx_norm")) + _bx_norm_name(isParamValid("bx_norm") ? getParam("bx_norm") : "") { #ifdef LIBMESH_HAVE_SLEPC if (_nl_sys_names.size() > 1) @@ -418,7 +418,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 / (_bx_norm_provided ? formNorm() : _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); @@ -642,8 +642,8 @@ EigenProblem::initPetscOutput() Real EigenProblem::formNorm() { - mooseAssert(_bx_norm_provided, + mooseAssert(!_bx_norm_name.empty(), "We should not get here unless a bx_norm postprocessor has been provided"); - return getPostprocessorValueByName(getParam("bx_norm")); + return getPostprocessorValueByName(_bx_norm_name); } #endif From 157c898b7ff8e5c7229dcce91129ba28ecdd62f2 Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Wed, 8 Nov 2023 15:17:02 -0700 Subject: [PATCH 11/16] Always respect required number of free power its --- framework/src/utils/SlepcSupport.C | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/framework/src/utils/SlepcSupport.C b/framework/src/utils/SlepcSupport.C index 9e2122e61be9..11b8ad922606 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 @@ -1103,8 +1107,7 @@ 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 From d892443d9eff2fabad4715c829a4a336e42ca870 Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Wed, 8 Nov 2023 15:17:33 -0700 Subject: [PATCH 12/16] Remove no longer needed tolerance specifications --- test/tests/executioners/eigen_convergence/a.i | 3 --- test/tests/executioners/eigen_convergence/tests | 5 ++--- 2 files changed, 2 insertions(+), 6 deletions(-) diff --git a/test/tests/executioners/eigen_convergence/a.i b/test/tests/executioners/eigen_convergence/a.i index f67516a5feb3..69a1aad33486 100644 --- a/test/tests/executioners/eigen_convergence/a.i +++ b/test/tests/executioners/eigen_convergence/a.i @@ -77,9 +77,6 @@ solve_type = PJFNK free_power_iterations = 4 nl_abs_tol = 2e-10 - # The first set of power iterations don't seem to respect these options - # petsc_options_iname = '-eps_tol' - # petsc_options_value = '1e-10' [] [Outputs] diff --git a/test/tests/executioners/eigen_convergence/tests b/test/tests/executioners/eigen_convergence/tests index 8d0b78d062ef..1848c1dd43e0 100644 --- a/test/tests/executioners/eigen_convergence/tests +++ b/test/tests/executioners/eigen_convergence/tests @@ -16,7 +16,6 @@ input = 'a.i' expect_out = '1[01] Nonlinear' absent_out = '12 Nonlinear' - cli_args = '-eps_tol 2e-10' 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] @@ -26,7 +25,7 @@ input = 'a.i' expect_out = '7 Nonlinear' absent_out = '8 Nonlinear' - cli_args = 'Problem/bx_norm=fluxintegral -eps_tol 2e-10' + 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.' [] @@ -37,7 +36,7 @@ input = 'a.i' expect_out = '4 Nonlinear' absent_out = '5 Nonlinear' - cli_args = 'Problem/bx_norm=fluxintegral -eps_power_sign_normalization 0 -eps_tol 2e-10' + 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.' [] From 058a85d4aca01b73bc6afef59ade67d6a7b21b40 Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Thu, 9 Nov 2023 10:19:59 -0800 Subject: [PATCH 13/16] Use std::optional for _bx_norm_name --- framework/include/problems/EigenProblem.h | 4 ++-- framework/src/problems/EigenProblem.C | 8 +++++--- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/framework/include/problems/EigenProblem.h b/framework/include/problems/EigenProblem.h index d0774674cad4..bae7c1d7f04e 100644 --- a/framework/include/problems/EigenProblem.h +++ b/framework/include/problems/EigenProblem.h @@ -232,7 +232,7 @@ class EigenProblem : public FEProblemBase /** * Whether a Bx norm postprocessor has been provided */ - bool bxNormProvided() const { return !_bx_norm_name.empty(); } + bool bxNormProvided() const { return _bx_norm_name.has_value(); } protected: unsigned int _n_eigen_pairs_required; @@ -284,7 +284,7 @@ class EigenProblem : public FEProblemBase /// 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 PostprocessorName _bx_norm_name; + const std::optional _bx_norm_name; #endif diff --git a/framework/src/problems/EigenProblem.C b/framework/src/problems/EigenProblem.C index 3348bf97a40c..cfce0f6d9d81 100644 --- a/framework/src/problems/EigenProblem.C +++ b/framework/src/problems/EigenProblem.C @@ -62,7 +62,9 @@ EigenProblem::EigenProblem(const InputParameters & parameters) _has_normalization(false), _normal_factor(1.0), _first_solve(declareRestartableData("first_solve", true)), - _bx_norm_name(isParamValid("bx_norm") ? getParam("bx_norm") : "") + _bx_norm_name(isParamValid("bx_norm") + ? std::make_optional(getParam("bx_norm")) + : std::nullopt) { #ifdef LIBMESH_HAVE_SLEPC if (_nl_sys_names.size() > 1) @@ -642,8 +644,8 @@ EigenProblem::initPetscOutput() Real EigenProblem::formNorm() { - mooseAssert(!_bx_norm_name.empty(), + mooseAssert(_bx_norm_name, "We should not get here unless a bx_norm postprocessor has been provided"); - return getPostprocessorValueByName(_bx_norm_name); + return getPostprocessorValueByName(*_bx_norm_name); } #endif From 19bb90c9fad49428ecb4a6680f5e6529baaed41d Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Thu, 9 Nov 2023 10:23:50 -0800 Subject: [PATCH 14/16] Edit param descriptions Co-authored-by: Guillaume Giudicelli --- framework/src/kernels/MassEigenKernel.C | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/framework/src/kernels/MassEigenKernel.C b/framework/src/kernels/MassEigenKernel.C index 0151a31c9e2c..dd5006040eaa 100644 --- a/framework/src/kernels/MassEigenKernel.C +++ b/framework/src/kernels/MassEigenKernel.C @@ -15,9 +15,10 @@ 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.addParam("coefficient", 1.0, "Coefficient of the term"); + 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; } From 16ded97a1c9d6e15b8b4191ed6a1ff6a09777fd9 Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Thu, 9 Nov 2023 10:47:17 -0800 Subject: [PATCH 15/16] Add comments to attachCallbacksToMat --- framework/src/utils/SlepcSupport.C | 23 ++++++++++++++++++++++- 1 file changed, 22 insertions(+), 1 deletion(-) diff --git a/framework/src/utils/SlepcSupport.C b/framework/src/utils/SlepcSupport.C index 11b8ad922606..2360c75ada56 100644 --- a/framework/src/utils/SlepcSupport.C +++ b/framework/src/utils/SlepcSupport.C @@ -978,21 +978,41 @@ mooseSlepcEigenFormNorm(SNES /*snes*/, Vec /*Bx*/, PetscReal * norm, void * ctx) 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); @@ -1107,7 +1127,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 From e99398534754d440700b6a134b9f460d51bed856 Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Thu, 9 Nov 2023 10:47:59 -0800 Subject: [PATCH 16/16] Remove what I believe should be unnecessary PetscObjectCompose calls --- framework/src/utils/SlepcSupport.C | 5 ----- 1 file changed, 5 deletions(-) diff --git a/framework/src/utils/SlepcSupport.C b/framework/src/utils/SlepcSupport.C index 2360c75ada56..a9d54b0e40ac 100644 --- a/framework/src/utils/SlepcSupport.C +++ b/framework/src/utils/SlepcSupport.C @@ -1016,15 +1016,10 @@ attachCallbacksToMat(EigenProblem & eigen_problem, Mat mat, bool eigen) 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", nullptr); PetscObjectCompose((PetscObject)mat, "formNormCtx", (PetscObject)container); - } PetscContainerDestroy(&container); }