From 98f76b013b1a85d44ea547c0db15bdb6764eef30 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Thu, 4 May 2023 10:21:15 -0600 Subject: [PATCH] Move functions into local lambdas, so that a member variable can be made non-public. --- include/deal.II/sundials/kinsol.h | 15 +- source/sundials/kinsol.cc | 343 ++++++++++++++---------------- 2 files changed, 170 insertions(+), 188 deletions(-) diff --git a/include/deal.II/sundials/kinsol.h b/include/deal.II/sundials/kinsol.h index 024ba92e22..4ac2c60247 100644 --- a/include/deal.II/sundials/kinsol.h +++ b/include/deal.II/sundials/kinsol.h @@ -699,14 +699,6 @@ namespace SUNDIALS << "returned a negative error code: " << arg1 << ". Please consult SUNDIALS manual."); - - /** - * A pointer to any exception that may have been thrown in user-defined - * call-backs and that we have to deal after the KINSOL function we call - * has returned. - */ - mutable std::exception_ptr pending_exception; - private: /** * Throw an exception when a function with the given name is not @@ -766,6 +758,13 @@ namespace SUNDIALS * Memory pool of vectors. */ GrowingVectorMemory mem; + + /** + * A pointer to any exception that may have been thrown in user-defined + * call-backs and that we have to deal after the KINSOL function we call + * has returned. + */ + mutable std::exception_ptr pending_exception; }; } // namespace SUNDIALS diff --git a/source/sundials/kinsol.cc b/source/sundials/kinsol.cc index dc57fc7a23..bdf44acd71 100644 --- a/source/sundials/kinsol.cc +++ b/source/sundials/kinsol.cc @@ -199,180 +199,6 @@ namespace SUNDIALS } - namespace - { - template - int - residual_callback(N_Vector yy, N_Vector FF, void *user_data) - { - KINSOL &solver = - *static_cast *>(user_data); - GrowingVectorMemory mem; - - typename VectorMemory::Pointer src_yy(mem); - solver.reinit_vector(*src_yy); - - typename VectorMemory::Pointer dst_FF(mem); - solver.reinit_vector(*dst_FF); - - internal::copy(*src_yy, yy); - - Assert(solver.residual, ExcInternalError()); - - const int err = call_and_possibly_capture_exception( - solver.residual, solver.pending_exception, *src_yy, *dst_FF); - - internal::copy(FF, *dst_FF); - - return err; - } - - - - template - int - iteration_callback(N_Vector yy, N_Vector FF, void *user_data) - { - KINSOL &solver = - *static_cast *>(user_data); - GrowingVectorMemory mem; - - typename VectorMemory::Pointer src_yy(mem); - solver.reinit_vector(*src_yy); - - typename VectorMemory::Pointer dst_FF(mem); - solver.reinit_vector(*dst_FF); - - internal::copy(*src_yy, yy); - - Assert(solver.iteration_function, ExcInternalError()); - - const int err = call_and_possibly_capture_exception( - solver.iteration_function, solver.pending_exception, *src_yy, *dst_FF); - - internal::copy(FF, *dst_FF); - - return err; - } - - - - template - int - setup_jacobian_callback(N_Vector u, - N_Vector f, - SUNMatrix /* ignored */, - void *user_data, - N_Vector /* tmp1 */, - N_Vector /* tmp2 */) - { - // Receive the object that describes the linear solver and - // unpack the pointer to the KINSOL object from which we can then - // get the 'setup' function. - const KINSOL &solver = - *static_cast *>(user_data); - - // Allocate temporary (deal.II-type) vectors into which to copy the - // N_vectors - GrowingVectorMemory mem; - typename VectorMemory::Pointer ycur(mem); - typename VectorMemory::Pointer fcur(mem); - solver.reinit_vector(*ycur); - solver.reinit_vector(*fcur); - - internal::copy(*ycur, u); - internal::copy(*fcur, f); - - // Call the user-provided setup function with these arguments: - return call_and_possibly_capture_exception(solver.setup_jacobian, - solver.pending_exception, - *ycur, - *fcur); - } - - - - template - int - solve_with_jacobian_callback(SUNLinearSolver LS, - SUNMatrix /*ignored*/, - N_Vector x, - N_Vector b, - realtype tol) - { - // Receive the object that describes the linear solver and - // unpack the pointer to the KINSOL object from which we can then - // get the 'reinit' and 'solve' functions. - const KINSOL &solver = - *static_cast *>(LS->content); - - // This is where we have to make a decision about which of the two - // signals to call. Let's first check the more modern one: - if (solver.solve_with_jacobian) - { - // Allocate temporary (deal.II-type) vectors into which to copy the - // N_vectors - GrowingVectorMemory mem; - typename VectorMemory::Pointer src_b(mem); - typename VectorMemory::Pointer dst_x(mem); - - solver.reinit_vector(*src_b); - solver.reinit_vector(*dst_x); - - internal::copy(*src_b, b); - - const int err = - call_and_possibly_capture_exception(solver.solve_with_jacobian, - solver.pending_exception, - *src_b, - *dst_x, - tol); - - internal::copy(x, *dst_x); - - return err; - } - else - { - // User has not provided the modern callback, so the fact that we are - // here means that they must have given us something for the old - // signal. Check this. - Assert(solver.solve_jacobian_system, ExcInternalError()); - - // Allocate temporary (deal.II-type) vectors into which to copy the - // N_vectors - GrowingVectorMemory mem; - typename VectorMemory::Pointer src_ycur(mem); - typename VectorMemory::Pointer src_fcur(mem); - typename VectorMemory::Pointer src_b(mem); - typename VectorMemory::Pointer dst_x(mem); - - solver.reinit_vector(*src_b); - solver.reinit_vector(*dst_x); - - internal::copy(*src_b, b); - - // Call the user-provided setup function with these arguments. Note - // that Sundials 4.x and later no longer provide values for - // src_ycur and src_fcur, and so we simply pass dummy vector in. - // These vectors will have zero lengths because we don't reinit them - // above. - const int err = - call_and_possibly_capture_exception(solver.solve_jacobian_system, - solver.pending_exception, - *src_ycur, - *src_fcur, - *src_b, - *dst_x); - - internal::copy(x, *dst_x); - - return err; - } - } - } // namespace - - template KINSOL::KINSOL(const AdditionalData &data) @@ -384,8 +210,7 @@ namespace SUNDIALS template KINSOL::KINSOL(const AdditionalData &data, const MPI_Comm & mpi_comm) - : pending_exception(nullptr) - , data(data) + : data(data) , mpi_communicator(mpi_comm) , kinsol_mem(nullptr) # if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) @@ -394,6 +219,7 @@ namespace SUNDIALS , solution(nullptr) , u_scale(nullptr) , f_scale(nullptr) + , pending_exception(nullptr) { set_functions_to_trigger_an_assert(); @@ -555,9 +381,62 @@ namespace SUNDIALS AssertKINSOL(status); if (data.strategy == AdditionalData::fixed_point) - status = KINInit(kinsol_mem, iteration_callback, solution); + status = KINInit( + kinsol_mem, + /* wrap up the iteration_function() callback: */ + [](N_Vector yy, N_Vector FF, void *user_data) -> int { + KINSOL &solver = + *static_cast *>(user_data); + GrowingVectorMemory mem; + + typename VectorMemory::Pointer src_yy(mem); + solver.reinit_vector(*src_yy); + + typename VectorMemory::Pointer dst_FF(mem); + solver.reinit_vector(*dst_FF); + + internal::copy(*src_yy, yy); + + Assert(solver.iteration_function, ExcInternalError()); + + const int err = + call_and_possibly_capture_exception(solver.iteration_function, + solver.pending_exception, + *src_yy, + *dst_FF); + + internal::copy(FF, *dst_FF); + + return err; + }, + solution); else - status = KINInit(kinsol_mem, residual_callback, solution); + status = KINInit( + kinsol_mem, + /* wrap up the residual() callback: */ + [](N_Vector yy, N_Vector FF, void *user_data) -> int { + KINSOL &solver = + *static_cast *>(user_data); + GrowingVectorMemory mem; + + typename VectorMemory::Pointer src_yy(mem); + solver.reinit_vector(*src_yy); + + typename VectorMemory::Pointer dst_FF(mem); + solver.reinit_vector(*dst_FF); + + internal::copy(*src_yy, yy); + + Assert(solver.residual, ExcInternalError()); + + const int err = call_and_possibly_capture_exception( + solver.residual, solver.pending_exception, *src_yy, *dst_FF); + + internal::copy(FF, *dst_FF); + + return err; + }, + solution); AssertKINSOL(status); status = KINSetFuncNormTol(kinsol_mem, data.function_tolerance); @@ -620,7 +499,81 @@ namespace SUNDIALS return 0; }; - LS->ops->solve = solve_with_jacobian_callback; + LS->ops->solve = [](SUNLinearSolver LS, + SUNMatrix /*ignored*/, + N_Vector x, + N_Vector b, + realtype tol) -> int { + // Receive the object that describes the linear solver and + // unpack the pointer to the KINSOL object from which we can then + // get the 'reinit' and 'solve' functions. + const KINSOL &solver = + *static_cast *>(LS->content); + + // This is where we have to make a decision about which of the two + // signals to call. Let's first check the more modern one: + if (solver.solve_with_jacobian) + { + // Allocate temporary (deal.II-type) vectors into which to copy + // the N_vectors + GrowingVectorMemory mem; + typename VectorMemory::Pointer src_b(mem); + typename VectorMemory::Pointer dst_x(mem); + + solver.reinit_vector(*src_b); + solver.reinit_vector(*dst_x); + + internal::copy(*src_b, b); + + const int err = + call_and_possibly_capture_exception(solver.solve_with_jacobian, + solver.pending_exception, + *src_b, + *dst_x, + tol); + + internal::copy(x, *dst_x); + + return err; + } + else + { + // User has not provided the modern callback, so the fact that + // we are here means that they must have given us something for + // the old signal. Check this. + Assert(solver.solve_jacobian_system, ExcInternalError()); + + // Allocate temporary (deal.II-type) vectors into which to copy + // the N_vectors + GrowingVectorMemory mem; + typename VectorMemory::Pointer src_ycur(mem); + typename VectorMemory::Pointer src_fcur(mem); + typename VectorMemory::Pointer src_b(mem); + typename VectorMemory::Pointer dst_x(mem); + + solver.reinit_vector(*src_b); + solver.reinit_vector(*dst_x); + + internal::copy(*src_b, b); + + // Call the user-provided setup function with these arguments. + // Note that Sundials 4.x and later no longer provide values for + // src_ycur and src_fcur, and so we simply pass dummy vector in. + // These vectors will have zero lengths because we don't reinit + // them above. + const int err = call_and_possibly_capture_exception( + solver.solve_jacobian_system, + solver.pending_exception, + *src_ycur, + *src_fcur, + *src_b, + *dst_x); + + internal::copy(x, *dst_x); + + return err; + } + }; // Even though we don't use it, KINSOL still wants us to set some // kind of matrix object for the nonlinear solver. This is because @@ -663,7 +616,37 @@ namespace SUNDIALS setup_jacobian = [](const VectorType &, const VectorType &) { return 0; }; - status = KINSetJacFn(kinsol_mem, &setup_jacobian_callback); + status = KINSetJacFn( + kinsol_mem, + [](N_Vector u, + N_Vector f, + SUNMatrix /* ignored */, + void *user_data, + N_Vector /* tmp1 */, + N_Vector /* tmp2 */) { + // Receive the object that describes the linear solver and + // unpack the pointer to the KINSOL object from which we can then + // get the 'setup' function. + const KINSOL &solver = + *static_cast *>(user_data); + + // Allocate temporary (deal.II-type) vectors into which to copy the + // N_vectors + GrowingVectorMemory mem; + typename VectorMemory::Pointer ycur(mem); + typename VectorMemory::Pointer fcur(mem); + solver.reinit_vector(*ycur); + solver.reinit_vector(*fcur); + + internal::copy(*ycur, u); + internal::copy(*fcur, f); + + // Call the user-provided setup function with these arguments: + return call_and_possibly_capture_exception(solver.setup_jacobian, + solver.pending_exception, + *ycur, + *fcur); + }); AssertKINSOL(status); } -- 2.39.5