diff --git a/src/poly/solver/private/polySolver.cpp b/src/poly/solver/private/polySolver.cpp index 767618c..60325da 100644 --- a/src/poly/solver/private/polySolver.cpp +++ b/src/poly/solver/private/polySolver.cpp @@ -114,24 +114,25 @@ void PolySolver::assembleBlockSystem() { const std::unique_ptr forms = buildIndividualForms(blockOffsets); - const double penalty_param = m_config.get("Poly::Solver::ZeroDerivativePenalty", 1.0); - mfem::Array thetaCenterDofs, phiCenterDofs; - std::tie(thetaCenterDofs, phiCenterDofs) = findCenterElement(); - mfem::SparseMatrix& D_mat = forms->D->SpMat(); - - for (int i = 0; i < phiCenterDofs.Size(); ++i) - { - const int dof_idx = phiCenterDofs[i]; - if (dof_idx >= 0 && dof_idx < D_mat.Height()) { - D_mat(dof_idx, dof_idx) += penalty_param; - } - } + // const double penalty_param = m_config.get("Poly::Solver::ZeroDerivativePenalty", 1.0); + // mfem::Array thetaCenterDofs, phiCenterDofs; + // std::tie(thetaCenterDofs, phiCenterDofs) = findCenterElement(); + // mfem::SparseMatrix& D_mat = forms->D->SpMat(); + // + // for (int i = 0; i < phiCenterDofs.Size(); ++i) + // { + // const int dof_idx = phiCenterDofs[i]; + // if (dof_idx >= 0 && dof_idx < D_mat.Height()) { + // D_mat(dof_idx, dof_idx) += penalty_param; + // } + // } // --- Build the BlockOperator --- m_polytropOperator = std::make_unique( std::move(forms->M), std::move(forms->Q), std::move(forms->D), + std::move(forms->S), std::move(forms->f), blockOffsets); } @@ -152,26 +153,26 @@ std::unique_ptr PolySolver::buildIndividualForms(const mfem::Array(m_fePhi.get(), m_feTheta.get()), std::make_unique(m_feTheta.get(), m_fePhi.get()), std::make_unique(m_fePhi.get()), + std::make_unique(m_feTheta.get()), std::make_unique(m_feTheta.get()) ); - // --- -M negation -> M --- - mfem::Vector negOneVec(m_mesh.SpaceDimension()); - negOneVec = -1.0; - - m_negationCoeff = std::make_unique(negOneVec); - // --- Add the integrators to the forms --- - forms->M->AddDomainIntegrator(new mfem::MixedVectorWeakDivergenceIntegrator()); - forms->Q->AddDomainIntegrator(new mfem::MixedVectorGradientIntegrator()); - forms->D->AddDomainIntegrator(new mfem::VectorFEMassIntegrator()); + forms->M->AddDomainIntegrator(new mfem::MixedVectorWeakDivergenceIntegrator()); // M ∫∇ψ^θ·N^φ dV + forms->Q->AddDomainIntegrator(new mfem::MixedVectorGradientIntegrator()); // Q ∫ψ^φ·∇N^θ dV + forms->D->AddDomainIntegrator(new mfem::VectorFEMassIntegrator()); // D ∫ψ^φ·N^φ dV + forms->S->AddDomainIntegrator(new mfem::DiffusionIntegrator()); // S ∫∇ψ^θ·∇N^θ dV + forms->f->AddDomainIntegrator(new polyMFEMUtils::NonlinearPowerIntegrator(m_polytropicIndex)); // --- Assemble and Finalize the forms --- assembleAndFinalizeForm(forms->M); assembleAndFinalizeForm(forms->Q); assembleAndFinalizeForm(forms->D); + assembleAndFinalizeForm(forms->S); + + // Note: The NonlinearForm does not need to / cannot be finalized, as it is not a matrix form. Rather, the operator + // will evaluate the nonlinear form during the solve phase. - forms->f->AddDomainIntegrator(new polyMFEMUtils::NonlinearPowerIntegrator(m_polytropicIndex)); return forms; } @@ -179,7 +180,7 @@ std::unique_ptr PolySolver::buildIndividualForms(const mfem::Array Computes the local element matrices across the domain. Adds these to the global matrix . Adds these to the global matrix. - // Finalize => Builds the sparsity pattern and allows the SparseMatrix representation to be extracted. + // Finalize => Builds the sparsity pattern and allows the SparseMatrix representation to be extracted (CSR). f->Assemble(); f->Finalize(); } diff --git a/src/poly/solver/public/polySolver.h b/src/poly/solver/public/polySolver.h index 8b14be6..3e140fd 100644 --- a/src/poly/solver/public/polySolver.h +++ b/src/poly/solver/public/polySolver.h @@ -49,6 +49,7 @@ struct formBundle { std::unique_ptr M; //<-- M ∫∇ψ^θ·N^φ dV std::unique_ptr Q; //<-- Q ∫ψ^φ·∇N^θ dV std::unique_ptr D; //<-- D ∫ψ^φ·N^φ dV + std::unique_ptr S; //<-- S ∫∇ψ^θ·∇N^θ dV std::unique_ptr f; //<-- f(θ) ∫ψ^θ·θ^n dV }; @@ -84,8 +85,6 @@ private: // Private Attributes std::unique_ptr m_prec; - std::unique_ptr m_negationCoeff; - private: // Private methods PolySolver(mfem::Mesh& mesh, double n, double order); diff --git a/src/poly/utils/private/polytropeOperator.cpp b/src/poly/utils/private/polytropeOperator.cpp index 898c03c..4b4ee82 100644 --- a/src/poly/utils/private/polytropeOperator.cpp +++ b/src/poly/utils/private/polytropeOperator.cpp @@ -134,10 +134,10 @@ PolytropeOperator::PolytropeOperator( std::unique_ptr M, std::unique_ptr Q, std::unique_ptr D, + std::unique_ptr S, std::unique_ptr f, const mfem::Array &blockOffsets) : - // TODO: Need to update this so that the size is that of the reduced system operator mfem::Operator(blockOffsets.Last()), // Initialize the base class with the total size of the block offset vector m_blockOffsets(blockOffsets), m_jacobian(nullptr) { @@ -145,6 +145,7 @@ PolytropeOperator::PolytropeOperator( m_M = std::move(M); m_Q = std::move(Q); m_D = std::move(D); + m_S = std::move(S); m_f = std::move(f); // Use Gauss-Seidel smoother to approximate the inverse of the matrix @@ -181,21 +182,33 @@ void PolytropeOperator::Mult(const mfem::Vector &xFree, mfem::Vector &yFree) con mfem::Vector Mphi_term(theta_size); mfem::Vector Dphi_term(phi_size); mfem::Vector Qtheta_term(phi_size); + mfem::Vector Stheta_term(theta_size); + // Calculate R0 and R1 terms - // R0 = f(θ) + Mɸ - // R1 = Dɸ - Qθ + // R0 = f(θ) + (1+c)Mɸ + cSθ + // R1 = (1+c)Dɸ - (1+c)Qθ MFEM_ASSERT(m_f.get() != nullptr, "NonlinearForm m_f is null in PolytropeOperator::Mult"); - m_f->Mult(x_theta, f_term); - m_M->Mult(x_phi, Mphi_term); - m_D->Mult(x_phi, Dphi_term); - m_Q->Mult(x_theta, Qtheta_term); + m_f->Mult(x_theta, f_term); // f(θ) + m_M->Mult(x_phi, Mphi_term); // Mɸ + m_D->Mult(x_phi, Dphi_term); // Dɸ + m_Q->Mult(x_theta, Qtheta_term); // Qθ + m_S->Mult(x_theta, Stheta_term); // Sθ - add(f_term, Mphi_term, y_R0); - subtract(Dphi_term, Qtheta_term, y_R1); + // PERF: these multiplications may be able to be refactored into the matrices so they only need to be done once. + Stheta_term *= m_stabilizationCoefficient; // cSθ + Qtheta_term *= m_IncrementedStabilizationCoefficient; // (1+c)Qθ + Mphi_term *= m_IncrementedStabilizationCoefficient; // (1+c)Mɸ + Dphi_term *= m_IncrementedStabilizationCoefficient; // (1+c)Dɸ + mfem::Vector y_R0_temp(theta_size); + add(f_term, Mphi_term, y_R0_temp); // R0 = f(θ) + (1+c)Mɸ + add(y_R0_temp, Stheta_term, y_R0); // R0 = f(θ) + (1+c)Mɸ + cSθ + subtract(Dphi_term, Qtheta_term, y_R1); // R1 = (1+c)Dɸ - (1+c)Qθ + + // --- Scatter the residual vector y into the free dofs --- yFree.SetSize(m_reducedBlockOffsets.Last()); MFEM_ASSERT(m_freeDofs.Size() == m_reducedBlockOffsets.Last(), "PolytropeOperator::Mult: Size of free dofs does not match reduced block offsets size."); for (int i = 0, j = 0; i < y.Size(); ++i) { @@ -207,7 +220,6 @@ void PolytropeOperator::Mult(const mfem::Vector &xFree, mfem::Vector &yFree) con } mfem::Operator& PolytropeOperator::GetGradient(const mfem::Vector &xFree) const { - //TODO: This now needs to be updated to deal with the reduced system size if (!m_isFinalized) { MFEM_ABORT("PolytropeOperator::GetGradient called before finalize"); } @@ -219,7 +231,7 @@ mfem::Operator& PolytropeOperator::GetGradient(const mfem::Vector &xFree) const // PERF: There are a lot of copies and loops here, probably performance could be gained by flattering some of these. auto &grad = m_f->GetGradient(x_theta); // updatePreconditioner(grad); - const auto gradMatrix = dynamic_cast(&grad); + const auto gradMatrix = dynamic_cast(&grad); // ∂f(θ)/∂θ if (gradMatrix == nullptr) { MFEM_ABORT("PolytropeOperator::GetGradient: Gradient is not a SparseMatrix."); @@ -231,7 +243,11 @@ mfem::Operator& PolytropeOperator::GetGradient(const mfem::Vector &xFree) const m_theta_ess_tdofs.first, m_theta_ess_tdofs.first ) - ); + ); // ∂f(θ)/∂θ (now reduced to only the free DOFs) + + + // TODO: Confirm this actually does what I want it to do + *m_gradReduced += *m_ScaledSReduced; // ∂f(θ)/∂θ + cS (reduced to only the free DOFs) m_jacobian->SetBlock(0, 0, m_gradReduced.get()); @@ -346,29 +362,73 @@ void PolytropeOperator::populate_free_dof_array() { // PolytropeOperator: Private Helper Methods - Construction and Setup void PolytropeOperator::construct_matrix_representations() { + // --- Construct the sparse matrix representations of M, Q, D, and S --- m_Mmat = std::make_unique(m_M->SpMat()); m_Qmat = std::make_unique(m_Q->SpMat()); m_Dmat = std::make_unique(m_D->SpMat()); + m_Smat = std::make_unique(m_S->SpMat()); - m_MReduced = std::make_unique(serif::utilities::build_reduced_matrix(*m_Mmat, m_phi_ess_tdofs.first, m_theta_ess_tdofs.first)); - m_QReduced = std::make_unique(serif::utilities::build_reduced_matrix(*m_Qmat, m_theta_ess_tdofs.first, m_phi_ess_tdofs.first)); - m_DReduced = std::make_unique(serif::utilities::build_reduced_matrix(*m_Dmat, m_phi_ess_tdofs.first, m_phi_ess_tdofs.first)); + // --- Reduce the matrices to only the free DOFs --- + m_MReduced = std::make_unique( + serif::utilities::build_reduced_matrix( + *m_Mmat, + m_phi_ess_tdofs.first, + m_theta_ess_tdofs.first) + ); + m_QReduced = std::make_unique( + serif::utilities::build_reduced_matrix( + *m_Qmat, + m_theta_ess_tdofs.first, + m_phi_ess_tdofs.first) + ); + m_DReduced = std::make_unique( + serif::utilities::build_reduced_matrix( + *m_Dmat, + m_phi_ess_tdofs.first, + m_phi_ess_tdofs.first) + ); + m_SReduced = std::make_unique( + serif::utilities::build_reduced_matrix( + *m_Smat, + m_theta_ess_tdofs.first, + m_theta_ess_tdofs.first) + ); - m_negQ_mat = std::make_unique(m_QReduced.get(), -1.0); + // --- Scale the reduced matrices by the stabilization coefficients --- + mfem::SparseMatrix MScaledTemp(*m_MReduced); // Create a temporary copy of the M matrix for scaling + MScaledTemp *= m_IncrementedStabilizationCoefficient; // Scale the M matrix by the incremented stabilization coefficient + m_MScaledReduced = std::make_unique(MScaledTemp); // Store the scaled M matrix so that it persists + + mfem::SparseMatrix QScaledTemp(*m_QReduced); + QScaledTemp *= m_IncrementedStabilizationCoefficient; // Scale the Q matrix by the incremented stabilization coefficient + m_QScaledReduced = std::make_unique(QScaledTemp); // Store the scaled Q matrix so that it persists + + mfem::SparseMatrix DRScaledTemp(*m_DReduced); + DRScaledTemp *= m_IncrementedStabilizationCoefficient; // Scale the D matrix by the incremented stabilization coefficient + m_DScaledReduced = std::make_unique(DRScaledTemp); // Store the scaled D matrix so that it persists + + // Scale the S matrix by the stabilization coefficient (so that the allocations only need to be done once) + mfem::SparseMatrix SScaledTemp(*m_SReduced); // Create a temporary copy of the S matrix for scaling + SScaledTemp *= m_stabilizationCoefficient; // Scale the S matrix by the stabilization coefficient + m_ScaledSReduced = std::make_unique(SScaledTemp); // Store the scaled S matrix so that it persists + + // --- Create the scaled operator for -(1+c)Q --- + m_negQ_mat = std::make_unique(m_QScaledReduced.get(), -1.0); } void PolytropeOperator::construct_reduced_block_offsets() { m_reducedBlockOffsets.SetSize(3); m_reducedBlockOffsets[0] = 0; // R0 block (theta) m_reducedBlockOffsets[1] = m_MReduced->Height(); // R1 block (theta) - m_reducedBlockOffsets[2] = m_QReduced->Height() + m_reducedBlockOffsets[1]; // R2 block (phi) + m_reducedBlockOffsets[2] = m_QReduced->Height() + m_reducedBlockOffsets[1]; // R2 block (phi + theta) } void PolytropeOperator::construct_jacobian_constant_terms() { m_jacobian = std::make_unique(m_reducedBlockOffsets); - m_jacobian->SetBlock(0, 1, m_MReduced.get()); //<- M (constant) - m_jacobian->SetBlock(1, 0, m_negQ_mat.get()); //<- -Q (constant) - m_jacobian->SetBlock(1, 1, m_DReduced.get()); //<- D (constant) + + m_jacobian->SetBlock(0, 1, m_MScaledReduced.get()); //<- (1+c)M (constant, reduced to only free DOFs) + m_jacobian->SetBlock(1, 0, m_negQ_mat.get()); //<- -(1+c)Q (constant, reduced to only free DOFs) + m_jacobian->SetBlock(1, 1, m_DScaledReduced.get()); //<- (1+c)D (constant, reduced to only free DOFs) } void PolytropeOperator::scatter_boundary_conditions() { diff --git a/src/poly/utils/public/polytropeOperator.h b/src/poly/utils/public/polytropeOperator.h index db332df..828214b 100644 --- a/src/poly/utils/public/polytropeOperator.h +++ b/src/poly/utils/public/polytropeOperator.h @@ -152,6 +152,7 @@ public: std::unique_ptr M, std::unique_ptr Q, std::unique_ptr D, + std::unique_ptr S, std::unique_ptr f, const mfem::Array &blockOffsets); @@ -286,19 +287,33 @@ private: std::unique_ptr m_M; ///< Bilinear form M, coupling θ and φ. std::unique_ptr m_Q; ///< Bilinear form Q, coupling φ and θ. std::unique_ptr m_D; ///< Bilinear form D, acting on φ. + std::unique_ptr m_S; std::unique_ptr m_f; ///< Nonlinear form f, acting on θ. // --- Full Matrix Representations (owned, derived from forms) --- std::unique_ptr m_Mmat; ///< Sparse matrix representation of M. std::unique_ptr m_Qmat; ///< Sparse matrix representation of Q. std::unique_ptr m_Dmat; ///< Sparse matrix representation of D. + std::unique_ptr m_Smat; // --- Reduced Matrix Representations (owned, after eliminating essential DOFs) --- std::unique_ptr m_MReduced; ///< Reduced M matrix (free DOFs only). std::unique_ptr m_QReduced; ///< Reduced Q matrix (free DOFs only). std::unique_ptr m_DReduced; ///< Reduced D matrix (free DOFs only). + std::unique_ptr m_SReduced; ///< Reduced S matrix (free DOFs only, used for least squares stabilization). mutable std::unique_ptr m_gradReduced; ///< Reduced gradient operator (G) for the nonlinear part f(θ). + // --- Scaled Reduced Matrix Representations (owned, after eliminating essential DOFs and scaling by stabilization coefficients) --- + std::unique_ptr m_MScaledReduced; ///< Scaled M matrix (free DOFs only, scaled by stabilization coefficient). + std::unique_ptr m_QScaledReduced; ///< Scaled Q matrix (free DOFs only, scaled by stabilization coefficient). + std::unique_ptr m_DScaledReduced; ///< Scaled D matrix (free DOFs only, scaled by stabilization coefficient). + std::unique_ptr m_ScaledSReduced; ///< Scaled S matrix (free DOFs only, scaled by stabilization coefficient). + + + // --- Stabilization Coefficient --- (Perhapses this should be a user parameter...) + static constexpr double m_stabilizationCoefficient = -1.0; ///< Stabilization coefficient for the system, used to more tightly couple ∇θ and φ. + static constexpr double m_IncrementedStabilizationCoefficient = 1.0 + m_stabilizationCoefficient; ///< 1 + Stabilization coefficient for the system, used to more tightly couple ∇θ and φ. + // --- State Vectors and DOF Management --- mutable mfem::Vector m_state; ///< Full state vector [θ, φ]^T, including essential DOFs. mfem::Array m_freeDofs; ///< Array of indices for free (non-essential) DOFs.