From 535f33ecb569f50690050560bbb0234838560db1 Mon Sep 17 00:00:00 2001 From: Luca Bertagna Date: Fri, 8 Nov 2024 09:33:24 -0700 Subject: [PATCH 01/13] Fix compiler warnings --- tests/unit/evaluators/GatherDistributedParameter.cpp | 1 - tests/unit/evaluators/ScatterResidual.cpp | 1 - 2 files changed, 2 deletions(-) diff --git a/tests/unit/evaluators/GatherDistributedParameter.cpp b/tests/unit/evaluators/GatherDistributedParameter.cpp index d6c1f221ce..cd346a131c 100644 --- a/tests/unit/evaluators/GatherDistributedParameter.cpp +++ b/tests/unit/evaluators/GatherDistributedParameter.cpp @@ -47,7 +47,6 @@ TEUCHOS_UNIT_TEST(evaluator_unit_tester, gatherDistributedParametersHessianVec) { using EvalType = PHAL::AlbanyTraits::HessianVec; using Scalar = EvalType::ScalarT; - using vec_str_pairs = std::vector>; auto comm = Albany::getDefaultComm(); diff --git a/tests/unit/evaluators/ScatterResidual.cpp b/tests/unit/evaluators/ScatterResidual.cpp index 2651cb9f1f..d746f555bf 100644 --- a/tests/unit/evaluators/ScatterResidual.cpp +++ b/tests/unit/evaluators/ScatterResidual.cpp @@ -401,7 +401,6 @@ TEUCHOS_UNIT_TEST(evaluator_unit_tester, scatterResidualHessianVecTensorRank1) { using EvalType = PHAL::AlbanyTraits::HessianVec; using Scalar = EvalType::ScalarT; - using vec_str_pairs = std::vector>; constexpr auto ALL = Kokkos::ALL(); constexpr auto ADD = Albany::CombineMode::ADD; From 11980207e2adc773323d3648fbd4c06e30a3460c Mon Sep 17 00:00:00 2001 From: Luca Bertagna Date: Fri, 10 Jan 2025 15:53:20 -0700 Subject: [PATCH 02/13] Moved AdvDiff problem to corePDEs --- src/corePDEs/CMakeLists.txt | 5 ++ .../evaluators/PHAL_AdvDiffResid.cpp | 0 .../evaluators/PHAL_AdvDiffResid.hpp | 3 +- .../evaluators/PHAL_AdvDiffResid_Def.hpp | 55 ++++++++++++------- .../problems/Albany_AdvDiffProblem.cpp | 0 .../problems/Albany_AdvDiffProblem.hpp | 0 .../problems/Albany_CoreProblemFactory.cpp | 7 +++ src/demoPDEs/CMakeLists.txt | 5 -- .../problems/Albany_DemoProblemFactory.cpp | 7 --- .../AdvDiff/CMakeLists.txt | 4 +- .../{demoPDEs => corePDEs}/AdvDiff/input.yaml | 5 +- tests/corePDEs/CMakeLists.txt | 2 + tests/demoPDEs/CMakeLists.txt | 1 - 13 files changed, 52 insertions(+), 42 deletions(-) rename src/{demoPDEs => corePDEs}/evaluators/PHAL_AdvDiffResid.cpp (100%) rename src/{demoPDEs => corePDEs}/evaluators/PHAL_AdvDiffResid.hpp (95%) rename src/{demoPDEs => corePDEs}/evaluators/PHAL_AdvDiffResid_Def.hpp (70%) rename src/{demoPDEs => corePDEs}/problems/Albany_AdvDiffProblem.cpp (100%) rename src/{demoPDEs => corePDEs}/problems/Albany_AdvDiffProblem.hpp (100%) rename tests/{demoPDEs => corePDEs}/AdvDiff/CMakeLists.txt (93%) rename tests/{demoPDEs => corePDEs}/AdvDiff/input.yaml (98%) diff --git a/src/corePDEs/CMakeLists.txt b/src/corePDEs/CMakeLists.txt index fd95a597fe..68d88fd43b 100644 --- a/src/corePDEs/CMakeLists.txt +++ b/src/corePDEs/CMakeLists.txt @@ -6,20 +6,25 @@ set(SOURCES + problems/Albany_AdvDiffProblem.cpp problems/Albany_CoreProblemFactory.cpp problems/Albany_HeatProblem.cpp problems/Albany_PopulateMesh.cpp problems/Albany_SideLaplacianProblem.cpp + evaluators/PHAL_AdvDiffResid.cpp evaluators/PHAL_HeatEqResid.cpp evaluators/PHAL_SideLaplacianResidual.cpp evaluators/PHAL_ThermalConductivity.cpp ) set(HEADERS + problems/Albany_AdvDiffProblem.hpp problems/Albany_CoreProblemFactory.hpp problems/Albany_HeatProblem.hpp problems/Albany_PopulateMesh.hpp problems/Albany_SideLaplacianProblem.hpp + evaluators/PHAL_AdvDiffResid.hpp + evaluators/PHAL_AdvDiffResid_Def.hpp evaluators/PHAL_HeatEqResid.hpp evaluators/PHAL_HeatEqResid_Def.hpp evaluators/PHAL_SideLaplacianResidual.hpp diff --git a/src/demoPDEs/evaluators/PHAL_AdvDiffResid.cpp b/src/corePDEs/evaluators/PHAL_AdvDiffResid.cpp similarity index 100% rename from src/demoPDEs/evaluators/PHAL_AdvDiffResid.cpp rename to src/corePDEs/evaluators/PHAL_AdvDiffResid.cpp diff --git a/src/demoPDEs/evaluators/PHAL_AdvDiffResid.hpp b/src/corePDEs/evaluators/PHAL_AdvDiffResid.hpp similarity index 95% rename from src/demoPDEs/evaluators/PHAL_AdvDiffResid.hpp rename to src/corePDEs/evaluators/PHAL_AdvDiffResid.hpp index 91262144bd..fcb08bccec 100644 --- a/src/demoPDEs/evaluators/PHAL_AdvDiffResid.hpp +++ b/src/corePDEs/evaluators/PHAL_AdvDiffResid.hpp @@ -46,8 +46,7 @@ class AdvDiffResid : public PHX::EvaluatorWithBaseImpl, PHX::MDField UDot; double mu; //viscosity coefficient - double a; //advection coefficient - double b; //advection coefficient + Teuchos::Array beta; // Advection coefficients bool useAugForm; //use augmented form? int formType; //augmented form type diff --git a/src/demoPDEs/evaluators/PHAL_AdvDiffResid_Def.hpp b/src/corePDEs/evaluators/PHAL_AdvDiffResid_Def.hpp similarity index 70% rename from src/demoPDEs/evaluators/PHAL_AdvDiffResid_Def.hpp rename to src/corePDEs/evaluators/PHAL_AdvDiffResid_Def.hpp index 26ef9199fb..b7473ae681 100644 --- a/src/demoPDEs/evaluators/PHAL_AdvDiffResid_Def.hpp +++ b/src/corePDEs/evaluators/PHAL_AdvDiffResid_Def.hpp @@ -4,6 +4,7 @@ // in the file "license.txt" in the top-level Albany directory // //*****************************************************************// +#include "Albany_StringUtils.hpp" #include "Teuchos_TestForException.hpp" #include "Phalanx_DataLayout.hpp" @@ -59,8 +60,12 @@ AdvDiffResid(const Teuchos::ParameterList& p) : vecDim = dims[2]; - a = bf_list->get("Advection a", 1.0); - b = bf_list->get("Advection b", 1.0); + beta = bf_list->get>("Advection Coefficient"); + TEUCHOS_TEST_FOR_EXCEPTION (static_cast(beta.size())!=numDims, std::logic_error, + "Error! The input array for 'Advection Coefficient' should have the same extent as the mesh num dims.\n" + " - beta: [" << util::join(beta,",") << "]\n" + " - num dims: " << numDims << "\n"); + mu = bf_list->get("Viscosity mu", 0.1); useAugForm = bf_list->get("Use Augmented Form", false); formType = bf_list->get("Augmented Form Type", 1); @@ -73,7 +78,7 @@ AdvDiffResid(const Teuchos::ParameterList& p) : "Invalid Augmented Form Type: " << formType << "; valid options are 1 and 2"); - std::cout << "a, b, mu: " << a << ", " << b << ", " << mu << std::endl; + std::cout << "beta, mu: [" << util::join(beta,",") << "], " << ", " << mu << std::endl; std::cout << " vecDim = " << vecDim << std::endl; std::cout << " numDims = " << numDims << std::endl; std::cout << "Augmented Form? " << useAugForm << std::endl; @@ -108,15 +113,23 @@ evaluateFields(typename Traits::EvalData workset) for (std::size_t node=0; node < numNodes; ++node) { for (std::size_t i=0; i Date: Wed, 8 Jan 2025 15:22:49 -0700 Subject: [PATCH 03/13] Add method to retrieve dxdp from the discretization --- src/disc/Albany_AbstractDiscretization.hpp | 3 +++ src/disc/Albany_AbstractMeshFieldAccessor.hpp | 3 +++ src/disc/Albany_ExtrudedDiscretization.cpp | 8 +++++++ src/disc/Albany_ExtrudedDiscretization.hpp | 1 + .../omegah/Albany_OmegahDiscretization.hpp | 4 ++++ .../omegah/Albany_OmegahMeshFieldAccessor.hpp | 7 ++++++ .../stk/Albany_MultiSTKFieldContainer.cpp | 13 +++++++++-- .../stk/Albany_MultiSTKFieldContainer.hpp | 4 ++++ .../stk/Albany_OrdinarySTKFieldContainer.cpp | 22 ++++++++++++++++++- .../stk/Albany_OrdinarySTKFieldContainer.hpp | 4 ++++ src/disc/stk/Albany_STKDiscretization.cpp | 12 ++++++---- src/disc/stk/Albany_STKDiscretization.hpp | 3 +++ 12 files changed, 77 insertions(+), 7 deletions(-) diff --git a/src/disc/Albany_AbstractDiscretization.hpp b/src/disc/Albany_AbstractDiscretization.hpp index 3cb264ed2a..395de97547 100644 --- a/src/disc/Albany_AbstractDiscretization.hpp +++ b/src/disc/Albany_AbstractDiscretization.hpp @@ -297,6 +297,9 @@ class AbstractDiscretization virtual void getSolutionMV(Thyra_MultiVector& soln, bool overlapped = false) const = 0; + virtual void + getSolutionDxDp(Thyra_MultiVector& dxdp, bool overlapped = false) const = 0; + virtual void getField(Thyra_Vector& field_vector, const std::string& field_name) const = 0; virtual void diff --git a/src/disc/Albany_AbstractMeshFieldAccessor.hpp b/src/disc/Albany_AbstractMeshFieldAccessor.hpp index 7ce13d4914..34e15d6987 100644 --- a/src/disc/Albany_AbstractMeshFieldAccessor.hpp +++ b/src/disc/Albany_AbstractMeshFieldAccessor.hpp @@ -54,6 +54,9 @@ class AbstractMeshFieldAccessor virtual void fillSolnMultiVector (Thyra_MultiVector& soln, const dof_mgr_ptr_t& sol_dof_mgr, const bool overlapped) = 0; + virtual void fillSolnSensitivity (Thyra_MultiVector& dxdp, + const dof_mgr_ptr_t& sol_dof_mgr, + const bool overlapped) = 0; // Write to mesh methods virtual void saveVector (const Thyra_Vector& field_vector, diff --git a/src/disc/Albany_ExtrudedDiscretization.cpp b/src/disc/Albany_ExtrudedDiscretization.cpp index 333e0eb3a0..8f4dc632a0 100644 --- a/src/disc/Albany_ExtrudedDiscretization.cpp +++ b/src/disc/Albany_ExtrudedDiscretization.cpp @@ -241,6 +241,14 @@ ExtrudedDiscretization::getSolutionMV( // solutionFieldContainer->fillSolnMultiVector(result, getDOFManager(), overlapped); } +void +ExtrudedDiscretization::getSolutionDxDp( + Thyra_MultiVector& /* result */, + const bool /* overlapped */) const +{ + throw NotYetImplemented("ExtrudedDiscretization::getSolutionDxDp"); +} + /*****************************************************************/ /*** Private functions follow. These are just used in above code */ /*****************************************************************/ diff --git a/src/disc/Albany_ExtrudedDiscretization.hpp b/src/disc/Albany_ExtrudedDiscretization.hpp index 50bf974fc9..c2e9aa00b4 100644 --- a/src/disc/Albany_ExtrudedDiscretization.hpp +++ b/src/disc/Albany_ExtrudedDiscretization.hpp @@ -107,6 +107,7 @@ class ExtrudedDiscretization : public AbstractDiscretization Teuchos::RCP getSolutionField (const bool overlapped = false) const override; void getSolutionMV (Thyra_MultiVector& result, bool overlapped) const override; + void getSolutionDxDp (Thyra_MultiVector& result, bool overlapped) const override; void getField (Thyra_Vector& field_vector, const std::string& field_name) const override; void setField (const Thyra_Vector& field_vector, diff --git a/src/disc/omegah/Albany_OmegahDiscretization.hpp b/src/disc/omegah/Albany_OmegahDiscretization.hpp index 2afb813038..09a7aa96c1 100644 --- a/src/disc/omegah/Albany_OmegahDiscretization.hpp +++ b/src/disc/omegah/Albany_OmegahDiscretization.hpp @@ -132,6 +132,10 @@ class OmegahDiscretization : public AbstractDiscretization } void getSolutionMV (Thyra_MultiVector& solution, bool overlapped) const override; + void getSolutionDxDp (Thyra_MultiVector& /* result */, bool /* overlapped */) const override + { + TEUCHOS_TEST_FOR_EXCEPTION(true,NotYetImplemented,"OmegahDiscretization::getSolutionDxDp"); + } void getField( Thyra_Vector& field_vector, const std::string& field_name) const override; diff --git a/src/disc/omegah/Albany_OmegahMeshFieldAccessor.hpp b/src/disc/omegah/Albany_OmegahMeshFieldAccessor.hpp index 1264a9c89b..070444395a 100644 --- a/src/disc/omegah/Albany_OmegahMeshFieldAccessor.hpp +++ b/src/disc/omegah/Albany_OmegahMeshFieldAccessor.hpp @@ -41,6 +41,13 @@ class OmegahMeshFieldAccessor : public AbstractMeshFieldAccessor TEUCHOS_TEST_FOR_EXCEPTION(true,NotYetImplemented,"OmegahMeshFieldAccessor::fillSolnMultiVector"); } + void fillSolnSensitivity(Thyra_MultiVector& /* dxdp */, + const Teuchos::RCP& /* solution_dof_mgr */, + const bool /* overlapped */) override + { + TEUCHOS_TEST_FOR_EXCEPTION(true,NotYetImplemented,"OmegahMeshFieldAccessor::fillSolnSensitivity"); + } + // Write to mesh methods void saveVector (const Thyra_Vector& field_vector, const std::string& field_name, diff --git a/src/disc/stk/Albany_MultiSTKFieldContainer.cpp b/src/disc/stk/Albany_MultiSTKFieldContainer.cpp index 8ca56e0428..ad7242b858 100644 --- a/src/disc/stk/Albany_MultiSTKFieldContainer.cpp +++ b/src/disc/stk/Albany_MultiSTKFieldContainer.cpp @@ -10,6 +10,7 @@ #include "Albany_Macros.hpp" #include "Albany_MultiSTKFieldContainer.hpp" #include "Albany_STKFieldContainerHelper.hpp" +#include "Albany_DiscretizationUtils.hpp" #include "Albany_ThyraUtils.hpp" #include "Albany_GlobalLocalIndexer.hpp" #include "Albany_StringUtils.hpp" @@ -281,6 +282,14 @@ fillSolnMultiVector ( Thyra_MultiVector& solution, } } +void MultiSTKFieldContainer:: +fillSolnSensitivity(Thyra_MultiVector& /* dxdp */, + const Teuchos::RCP& /* solution_dof_mgr */, + const bool /* overlapped */) +{ + throw NotYetImplemented ("MultiSTKFieldContainer::fillSolnSensitivity"); +} + void MultiSTKFieldContainer:: saveVector (const Thyra_Vector& field_vector, const std::string& field_name, @@ -322,7 +331,7 @@ saveSolnVector (const Thyra_Vector& soln, { // TODO: why can't we save also solution_dot? auto out = Teuchos::VerboseObjectBase::getDefaultOStream(); - *out << "IKT WARNING: calling MultiSTKFieldContainer::saveSolnVectorT with soln_dot,\n" + *out << "IKT WARNING: calling MultiSTKFieldContainer::saveSolnVector with soln_dot,\n" << "but this function has not been extended to write soln_dot properly to the Exodus file.\n" << "Exodus file will contain only soln, not soln_dot.\n"; @@ -339,7 +348,7 @@ saveSolnVector (const Thyra_Vector& soln, { // TODO: why can't we save also solution_dot? auto out = Teuchos::VerboseObjectBase::getDefaultOStream(); - *out << "IKT WARNING: calling MultiSTKFieldContainer::saveSolnVectorT with soln_dot and soln_dotdot,\n" + *out << "IKT WARNING: calling MultiSTKFieldContainer::saveSolnVector with soln_dot and soln_dotdot,\n" << "but this function has not been extended to write soln_dot[dot] properly to the Exodus file.\n" << "Exodus file will contain only soln, not soln_dot nor soln_dotdot.\n"; diff --git a/src/disc/stk/Albany_MultiSTKFieldContainer.hpp b/src/disc/stk/Albany_MultiSTKFieldContainer.hpp index 06d0f88803..0300ecedd1 100644 --- a/src/disc/stk/Albany_MultiSTKFieldContainer.hpp +++ b/src/disc/stk/Albany_MultiSTKFieldContainer.hpp @@ -81,6 +81,10 @@ class MultiSTKFieldContainer : public GenericSTKFieldContainer const dof_mgr_ptr_t& sol_dof_mgr, const bool overlapped); + void fillSolnSensitivity (Thyra_MultiVector& dxdp, + const dof_mgr_ptr_t& soln_dof_mgr, + const bool overlapped); + void transferSolutionToCoords(); private: diff --git a/src/disc/stk/Albany_OrdinarySTKFieldContainer.cpp b/src/disc/stk/Albany_OrdinarySTKFieldContainer.cpp index c60d052166..5f759d1251 100644 --- a/src/disc/stk/Albany_OrdinarySTKFieldContainer.cpp +++ b/src/disc/stk/Albany_OrdinarySTKFieldContainer.cpp @@ -242,7 +242,27 @@ fillSolnMultiVector(Thyra_MultiVector& solution, TEUCHOS_TEST_FOR_EXCEPTION (col.is_null(), std::runtime_error, "Error! Could not extract column from multivector.\n"); - fillVectorImpl(*solution.col(icomp), solution_field[icomp]->name(), solution_dof_mgr, overlapped); + fillVectorImpl(*col, solution_field[icomp]->name(), solution_dof_mgr, overlapped); + } +} + +void OrdinarySTKFieldContainer:: +fillSolnSensitivity(Thyra_MultiVector& dxdp, + const Teuchos::RCP& solution_dof_mgr, + const bool overlapped) +{ + TEUCHOS_TEST_FOR_EXCEPTION( + dxdp.domain()->dim() != this->num_params, std::runtime_error, + "Error in fillSolnSensitivity! Wrong number of vectors in dxdp.\n" + " - num vectors: " << dxdp.domain()->dim() << "\n" + " - num_params : " << this->num_params << "\n"); + + for (int iparam = 0; iparam < solution_field_dxdp.size(); ++iparam) { + auto col = dxdp.col(iparam); + TEUCHOS_TEST_FOR_EXCEPTION (col.is_null(), std::runtime_error, + "Error! Could not extract column from multivector.\n"); + + fillVectorImpl(*col, solution_field_dxdp[iparam]->name(), solution_dof_mgr, overlapped); } } diff --git a/src/disc/stk/Albany_OrdinarySTKFieldContainer.hpp b/src/disc/stk/Albany_OrdinarySTKFieldContainer.hpp index c027d34e45..fed8852199 100644 --- a/src/disc/stk/Albany_OrdinarySTKFieldContainer.hpp +++ b/src/disc/stk/Albany_OrdinarySTKFieldContainer.hpp @@ -78,6 +78,10 @@ class OrdinarySTKFieldContainer : public GenericSTKFieldContainer const dof_mgr_ptr_t& soln_dof_mgr, const bool overlapped); + void fillSolnSensitivity (Thyra_MultiVector& dxdp, + const dof_mgr_ptr_t& soln_dof_mgr, + const bool overlapped); + void saveVector (const Thyra_Vector& field_vector, const std::string& field_name, const dof_mgr_ptr_t& field_dof_mgr, diff --git a/src/disc/stk/Albany_STKDiscretization.cpp b/src/disc/stk/Albany_STKDiscretization.cpp index 07b1c361b2..6fc4c4c746 100644 --- a/src/disc/stk/Albany_STKDiscretization.cpp +++ b/src/disc/stk/Albany_STKDiscretization.cpp @@ -996,8 +996,6 @@ STKDiscretization::getField(Thyra_Vector& result, const std::string& name) const void STKDiscretization::getSolutionField(Thyra_Vector& result, const bool overlapped) const { - TEUCHOS_TEST_FOR_EXCEPTION(overlapped, std::logic_error, "Not implemented."); - solutionFieldContainer->fillSolnVector(result, getDOFManager(), overlapped); } @@ -1006,11 +1004,17 @@ STKDiscretization::getSolutionMV( Thyra_MultiVector& result, const bool overlapped) const { - TEUCHOS_TEST_FOR_EXCEPTION(overlapped, std::logic_error, "Not implemented."); - solutionFieldContainer->fillSolnMultiVector(result, getDOFManager(), overlapped); } +void +STKDiscretization::getSolutionDxDp( + Thyra_MultiVector& result, + const bool overlapped) const +{ + solutionFieldContainer->fillSolnSensitivity(result, getDOFManager(), overlapped); +} + /*****************************************************************/ /*** Private functions follow. These are just used in above code */ /*****************************************************************/ diff --git a/src/disc/stk/Albany_STKDiscretization.hpp b/src/disc/stk/Albany_STKDiscretization.hpp index 20c349c456..eb2c7be1cb 100644 --- a/src/disc/stk/Albany_STKDiscretization.hpp +++ b/src/disc/stk/Albany_STKDiscretization.hpp @@ -197,6 +197,9 @@ class STKDiscretization : public AbstractDiscretization void getSolutionMV(Thyra_MultiVector& result, bool overlapped) const; + void + getSolutionDxDp (Thyra_MultiVector& result, bool overlapped) const; + void getField(Thyra_Vector& field_vector, const std::string& field_name) const; void From cabd5553e3c157fb51dc22dbb6cf616de2c986e0 Mon Sep 17 00:00:00 2001 From: Luca Bertagna Date: Wed, 8 Jan 2025 15:24:49 -0700 Subject: [PATCH 04/13] Fix nano-bug in STK disc --- src/disc/stk/Albany_STKDiscretization.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/disc/stk/Albany_STKDiscretization.cpp b/src/disc/stk/Albany_STKDiscretization.cpp index 6fc4c4c746..81259c9522 100644 --- a/src/disc/stk/Albany_STKDiscretization.cpp +++ b/src/disc/stk/Albany_STKDiscretization.cpp @@ -1077,7 +1077,8 @@ void STKDiscretization::computeVectorSpaces() it.second = create_dof_mgr(part_name, nodes_dof_name(), FE_Type::HGRAD,1,1); } - coordinates.resize(3 * getLocalSubdim(getOverlapNodeVectorSpace())); + const int meshDim = stkMeshStruct->numDim; + coordinates.resize(meshDim * getLocalSubdim(getOverlapNodeVectorSpace())); } void From a3cbe7b396ba5d9d737244a3cb68d04f5b6135c0 Mon Sep 17 00:00:00 2001 From: Luca Bertagna Date: Wed, 8 Jan 2025 15:28:11 -0700 Subject: [PATCH 05/13] Add disc interfaces for mesh adaptation --- src/disc/Albany_AbstractDiscretization.hpp | 10 +++++++++ src/disc/Albany_DiscretizationUtils.hpp | 22 ++++++++++++++++++- src/disc/Albany_ExtrudedDiscretization.cpp | 19 ++++++++++++++++ src/disc/Albany_ExtrudedDiscretization.hpp | 8 +++++++ .../omegah/Albany_OmegahDiscretization.hpp | 14 ++++++++++++ src/disc/stk/Albany_STKDiscretization.cpp | 16 ++++++++++++++ src/disc/stk/Albany_STKDiscretization.hpp | 20 ++++++++++++----- 7 files changed, 102 insertions(+), 7 deletions(-) diff --git a/src/disc/Albany_AbstractDiscretization.hpp b/src/disc/Albany_AbstractDiscretization.hpp index 395de97547..6ceb02d8f0 100644 --- a/src/disc/Albany_AbstractDiscretization.hpp +++ b/src/disc/Albany_AbstractDiscretization.hpp @@ -392,6 +392,16 @@ class AbstractDiscretization const bool overlapped = false, const bool force_write_solution = false) = 0; + // Check if mesh adaptation is needed, and if so what kind (topological or just mesh-movement) + virtual Teuchos::RCP + checkForAdaptation (const Teuchos::RCP& solution, + const Teuchos::RCP& solution_dot, + const Teuchos::RCP& solution_dotdot, + const Teuchos::RCP& dxdp) const = 0; + + // Check if mesh adaptation is needed, and if so adapt mesh (and possibly reinterpolate solution) + virtual void adapt (const Teuchos::RCP& adaptData) = 0; + protected: strmap_t> sideSetDiscretizations; diff --git a/src/disc/Albany_DiscretizationUtils.hpp b/src/disc/Albany_DiscretizationUtils.hpp index db3078f3ef..415d6a35eb 100644 --- a/src/disc/Albany_DiscretizationUtils.hpp +++ b/src/disc/Albany_DiscretizationUtils.hpp @@ -8,9 +8,10 @@ #define ALBANY_DISCRETIZATION_UTILS_HPP #include "Albany_KokkosTypes.hpp" +#include "Albany_ThyraTypes.hpp" #include "Albany_ScalarOrdinalTypes.hpp" -#include "Kokkos_DualView.hpp" +#include "Kokkos_DualView.hpp" #include "Intrepid2_Basis.hpp" #include "Teuchos_ArrayRCP.hpp" @@ -45,6 +46,25 @@ enum class MeshType { Unstructured // No structure known (e.g., read from file) }; +// When adapting mesh, we'll return an enum stating what kind of adaptation happened +enum class AdaptationType { + None, + Movement, + Topology +}; + +struct AdaptationData { + AdaptationType type = AdaptationType::None; + + // Current value of solution and its time deriv. + // If adapting, the discretization can interpolate these onto + // the new mesh. + Teuchos::RCP x; + Teuchos::RCP x_dot; + Teuchos::RCP x_dotdot; + Teuchos::RCP dxdp; +}; + inline std::string e2str (const FE_Type fe_type) { std::string s; diff --git a/src/disc/Albany_ExtrudedDiscretization.cpp b/src/disc/Albany_ExtrudedDiscretization.cpp index 8f4dc632a0..4b7b09b0fe 100644 --- a/src/disc/Albany_ExtrudedDiscretization.cpp +++ b/src/disc/Albany_ExtrudedDiscretization.cpp @@ -199,6 +199,25 @@ ExtrudedDiscretization::writeSolutionMVToFile( throw NotYetImplemented("ExtrudedDiscretization::writeSolutionMVToFile"); } +Teuchos::RCP +ExtrudedDiscretization:: +checkForAdaptation (const Teuchos::RCP& /* solution */, + const Teuchos::RCP& /* solution_dot */, + const Teuchos::RCP& /* solution_dotdot */, + const Teuchos::RCP& /* dxdp */) const +{ + // We can't just do + // return m_basal_disc->checkForAdaptation(solution,solution_dot,solution_dotdot); + // We need to decide WHAT to bass to basal disc: the whole solution or the projection? + throw NotYetImplemented("ExtrudedDiscretization::checkForAdaptation"); +} + +void ExtrudedDiscretization:: +adapt (const Teuchos::RCP& /* adaptData */) +{ + throw NotYetImplemented("ExtrudedDiscretization::adapt"); +} + Teuchos::RCP ExtrudedDiscretization::getSolutionField(bool /* overlapped */) const { diff --git a/src/disc/Albany_ExtrudedDiscretization.hpp b/src/disc/Albany_ExtrudedDiscretization.hpp index c2e9aa00b4..a392b00325 100644 --- a/src/disc/Albany_ExtrudedDiscretization.hpp +++ b/src/disc/Albany_ExtrudedDiscretization.hpp @@ -172,6 +172,14 @@ class ExtrudedDiscretization : public AbstractDiscretization const bool overlapped = false, const bool force_write_solution = false) override; + Teuchos::RCP + checkForAdaptation (const Teuchos::RCP& solution, + const Teuchos::RCP& solution_dot, + const Teuchos::RCP& solution_dotdot, + const Teuchos::RCP& dxdp) const override; + + void adapt (const Teuchos::RCP& adaptData) override; + void setFieldData(const Teuchos::RCP& sis) override; protected: diff --git a/src/disc/omegah/Albany_OmegahDiscretization.hpp b/src/disc/omegah/Albany_OmegahDiscretization.hpp index 09a7aa96c1..35999f1e64 100644 --- a/src/disc/omegah/Albany_OmegahDiscretization.hpp +++ b/src/disc/omegah/Albany_OmegahDiscretization.hpp @@ -243,6 +243,20 @@ class OmegahDiscretization : public AbstractDiscretization TEUCHOS_TEST_FOR_EXCEPTION(true,NotYetImplemented,"OmegahDiscretization::writeSolutionMVToFile"); } + Teuchos::RCP + checkForAdaptation (const Teuchos::RCP& /* solution */, + const Teuchos::RCP& /* solution_dot */, + const Teuchos::RCP& /* solution_dotdot */, + const Teuchos::RCP& /* dxdp */) const override + { + throw NotYetImplemented("OmegaDiscretization::checkForAdaptation"); + } + + void adapt (const Teuchos::RCP& /* adaptData */) override + { + throw NotYetImplemented("OmegaDiscretization::adapt"); + } + protected: Teuchos::RCP diff --git a/src/disc/stk/Albany_STKDiscretization.cpp b/src/disc/stk/Albany_STKDiscretization.cpp index 81259c9522..42741461f0 100644 --- a/src/disc/stk/Albany_STKDiscretization.cpp +++ b/src/disc/stk/Albany_STKDiscretization.cpp @@ -2442,4 +2442,20 @@ create_dof_mgr (const std::string& part_name, return dof_mgr; } +Teuchos::RCP +STKDiscretization:: +checkForAdaptation (const Teuchos::RCP& /* solution */, + const Teuchos::RCP& /* solution_dot */, + const Teuchos::RCP& /* solution_dotdot */, + const Teuchos::RCP& /* dxdp */) const +{ + return Teuchos::rcp(new AdaptationData()); +} + +void STKDiscretization:: +adapt (const Teuchos::RCP& /* adaptData */) +{ + throw NotYetImplemented ("STKDiscretization::adapt"); +} + } // namespace Albany diff --git a/src/disc/stk/Albany_STKDiscretization.hpp b/src/disc/stk/Albany_STKDiscretization.hpp index eb2c7be1cb..45f87d5b2c 100644 --- a/src/disc/stk/Albany_STKDiscretization.hpp +++ b/src/disc/stk/Albany_STKDiscretization.hpp @@ -285,13 +285,21 @@ class STKDiscretization : public AbstractDiscretization const bool overlapped = false, const bool force_write_solution = false); - /** Add a solution field - */ - void addSolutionField(const std::string & fieldName,const std::string & blockId); + Teuchos::RCP + checkForAdaptation (const Teuchos::RCP& solution, + const Teuchos::RCP& solution_dot, + const Teuchos::RCP& solution_dotdot, + const Teuchos::RCP& dxdp) const override; - /** Add a solution field - */ - void addCellField(const std::string & fieldName,const std::string & blockId); + void adapt (const Teuchos::RCP& adaptData) override; + + /** Add a solution field + */ + void addSolutionField(const std::string & fieldName,const std::string & blockId); + + /** Add a solution field + */ + void addCellField(const std::string & fieldName,const std::string & blockId); void setFieldData(const Teuchos::RCP& sis); From 7f84b400e97abb6fac4a0523ddcaf16b1273ee7c Mon Sep 17 00:00:00 2001 From: Luca Bertagna Date: Wed, 8 Jan 2025 15:29:03 -0700 Subject: [PATCH 06/13] Upgrade to SolutionManager * Store DxDp so it can be retrieved * Allow to reset vector spaces at runtime (for adaptation) --- src/SolutionManager.cpp | 64 +++++++++++++++++++++++++---------------- src/SolutionManager.hpp | 6 +++- 2 files changed, 44 insertions(+), 26 deletions(-) diff --git a/src/SolutionManager.cpp b/src/SolutionManager.cpp index 14a9841662..0d4a5a7b5b 100644 --- a/src/SolutionManager.cpp +++ b/src/SolutionManager.cpp @@ -34,7 +34,7 @@ SolutionManager::SolutionManager( Teuchos::RCP problemParams = Teuchos::sublist(appParams_, "Problem", true); - num_params_ = Albany::CalculateNumberParams(problemParams); + Albany::CalculateNumberParams(problemParams,&num_params_); // Want the initial time in the parameter library to be correct // if this is a restart solution @@ -66,31 +66,10 @@ SolutionManager::SolutionManager( paramLib_->setRealValue("Time", initialValue); } - { - auto owned_vs = disc_->getVectorSpace(); - auto overlapped_vs = disc_->getOverlapVectorSpace(); + reset_solution_space(true); - overlapped_soln = Thyra::createMembers(overlapped_vs, num_time_deriv + 1); - if (num_params_ > 0) { - overlapped_soln_dxdp = Thyra::createMembers(overlapped_vs, num_params_); - } else { - overlapped_soln_dxdp = Teuchos::null; - } - - // TODO: ditch the overlapped_*T and keep only overlapped_*. - // You need to figure out how to pass the graph in a Tpetra-free way - // though... - overlapped_f = Thyra::createMember(overlapped_vs); - - // This call allocates the non-overlapped MV - current_soln = Thyra::createMembers(owned_vs,num_time_deriv+1); - if (disc_->hasRestartSolution()) { - disc_->getSolutionMV(*current_soln); - } - - // Create the CombineAndScatterManager for handling distributed memory linear - // algebra communications - cas_manager = Albany::createCombineAndScatterManager(owned_vs, overlapped_vs); + if (disc_->hasRestartSolution()) { + disc_->getSolutionMV(*current_soln); } auto pbParams = Teuchos::sublist(appParams_, "Problem", true); @@ -143,6 +122,41 @@ SolutionManager::SolutionManager( } } +void SolutionManager:: +reset_solution_space (const bool initial_condition) +{ + auto owned_vs = disc_->getVectorSpace(); + auto overlapped_vs = disc_->getOverlapVectorSpace(); + + overlapped_soln = Thyra::createMembers(overlapped_vs, num_time_deriv + 1); + if (num_params_ > 0) { + overlapped_soln_dxdp = Thyra::createMembers(overlapped_vs, num_params_); + } else { + overlapped_soln_dxdp = Teuchos::null; + } + + // TODO: ditch the overlapped_*T and keep only overlapped_*. + // You need to figure out how to pass the graph in a Tpetra-free way + // though... + overlapped_f = Thyra::createMember(overlapped_vs); + + // This call allocates the non-overlapped MV + current_soln = Thyra::createMembers(owned_vs,num_time_deriv+1); + current_dxdp = Thyra::createMembers(owned_vs,num_params_); + + // Create the CombineAndScatterManager for handling distributed memory linear + // algebra communications + cas_manager = Albany::createCombineAndScatterManager(owned_vs, overlapped_vs); + + if (not initial_condition) { + // This must be a sol mgr reset after mesh adaptation. Get new solution from the disc + disc_->getSolutionMV(*current_soln); + if (num_params_ > 0) { + disc_->getSolutionDxDp(*current_dxdp); + } + } +} + Teuchos::RCP SolutionManager::updateAndReturnOverlapSolution( const Thyra_Vector& solution /* not overlapped */) diff --git a/src/SolutionManager.hpp b/src/SolutionManager.hpp index 3eb95b53e2..1927382dcb 100644 --- a/src/SolutionManager.hpp +++ b/src/SolutionManager.hpp @@ -42,7 +42,8 @@ class SolutionManager { Teuchos::RCP get_cas_manager() const { return cas_manager; } - Teuchos::RCP getCurrentSolution() { return current_soln; } + Teuchos::RCP getCurrentSolution() { return current_soln; } + Teuchos::RCP getCurrentDxDp() { return current_dxdp; } void scatterX (const Thyra_Vector& x, const Teuchos::Ptr x_dot, @@ -52,6 +53,7 @@ class SolutionManager { void scatterX (const Thyra_MultiVector& soln, const Teuchos::Ptr dxdp = Teuchos::null); + void reset_solution_space (const bool initial_condition); private: Teuchos::RCP cas_manager; @@ -60,6 +62,8 @@ class SolutionManager { // The solution directly from the discretization class Teuchos::RCP current_soln; + Teuchos::RCP current_dxdp; + Teuchos::RCP overlapped_soln; Teuchos::RCP overlapped_soln_dxdp; From 8299b5faa7d82b3154214710b07690a81feb18e6 Mon Sep 17 00:00:00 2001 From: Luca Bertagna Date: Wed, 8 Jan 2025 15:33:24 -0700 Subject: [PATCH 07/13] Add calls to adaptivity hooks for time-dep problems * Add new observer that is Tempus-compatible * Call adaptivity hooks in the new observer --- src/Albany_ModelEvaluator.cpp | 1 - src/Albany_PiroObserver.hpp | 2 +- src/Albany_PiroTempusObserver.cpp | 93 +++++++++++++++++++++++++++++++ src/Albany_PiroTempusObserver.hpp | 31 +++++++++++ src/Albany_SolverFactory.cpp | 9 ++- src/CMakeLists.txt | 1 + 6 files changed, 134 insertions(+), 3 deletions(-) create mode 100644 src/Albany_PiroTempusObserver.cpp create mode 100644 src/Albany_PiroTempusObserver.hpp diff --git a/src/Albany_ModelEvaluator.cpp b/src/Albany_ModelEvaluator.cpp index f674d753a8..48a69f1e2a 100644 --- a/src/Albany_ModelEvaluator.cpp +++ b/src/Albany_ModelEvaluator.cpp @@ -1113,7 +1113,6 @@ evalModelImpl(const Thyra_InArgs& inArgs, // Get the input arguments // - ObserverImpl observer(app); auto out = Teuchos::VerboseObjectBase::getDefaultOStream(); // Thyra vectors diff --git a/src/Albany_PiroObserver.hpp b/src/Albany_PiroObserver.hpp index 15f01adc71..1b599950f8 100644 --- a/src/Albany_PiroObserver.hpp +++ b/src/Albany_PiroObserver.hpp @@ -71,7 +71,7 @@ class PiroObserver : public Piro::ObserverBase { void parameterChanged (const std::string& param) { impl_.parameterChanged(param); } -private: +protected: void observeSolutionImpl (const Teuchos::RCP& x, const Teuchos::RCP& x_dot, diff --git a/src/Albany_PiroTempusObserver.cpp b/src/Albany_PiroTempusObserver.cpp new file mode 100644 index 0000000000..e8b6767b00 --- /dev/null +++ b/src/Albany_PiroTempusObserver.cpp @@ -0,0 +1,93 @@ +//*****************************************************************// +// Albany 3.0: Copyright 2016 Sandia Corporation // +// This Software is released under the BSD license detailed // +// in the file "license.txt" in the top-level Albany directory // +//*****************************************************************// + +#include "Albany_PiroTempusObserver.hpp" + +#include +#include +#include + +namespace Albany +{ + +PiroTempusObserver:: +PiroTempusObserver(const Teuchos::RCP& app, + const Teuchos::RCP& model) + : PiroObserver(app,model) + , app_ (app) +{ + // Nothing else to do +} + +void PiroTempusObserver:: +observeEndTimeStep(const Tempus::Integrator& integrator) +{ + auto& integrator_nc = const_cast&>(integrator); + auto history_nc = integrator_nc.getNonConstSolutionHistory(); + auto state_nc = history_nc->getCurrentState(); + + TEUCHOS_TEST_FOR_EXCEPTION (state_nc.is_null(), std::runtime_error, + "Error! Unexpectedly found a null current state in the tempus integrator.\n"); + + //Don't observe solution if step failed to converge + if (state_nc->getSolutionStatus() == Tempus::Status::FAILED) { + return; + } + + // In order for the DISC to decide whether to adapt or not, + // we need to write the solution in the mesh db. + // HOWEVER, we don't want to do a regular "observation" step, + // since we don't want to write the solution to file, or to observe responses + Teuchos::RCP dxdp; + auto integrator_ptr = Teuchos::rcpFromRef(integrator); + auto fwd_sens_integrator = Teuchos::rcp_dynamic_cast>(integrator_ptr); + if (Teuchos::nonnull(fwd_sens_integrator)) { + dxdp = fwd_sens_integrator->getDxDp(); + } + + auto time = state_nc->getTime(); + auto disc = app_->getDiscretization(); + auto adaptData = disc->checkForAdaptation(state_nc->getX(),state_nc->getXDot(),state_nc->getXDotDot(),dxdp); + + // Before observing the solution, check if we need to adapt + if (adaptData->type!=AdaptationType::None) { + disc->adapt (adaptData); + // Make the solution manager import the new solution from the discretization + app_->getAdaptSolMgr()->reset_solution_space(false); + auto num_time_derivs = app_->getNumTimeDerivs(); + + // Get new solution + auto sol = app_->getAdaptSolMgr()->getCurrentSolution(); + + // Reset vectors now that they have been adapted + state_nc->setX(sol->col(0)); + if (num_time_derivs>0) { + state_nc->setXDot(sol->col(1)); + if (num_time_derivs>1) { + state_nc->setXDotDot(sol->col(2)); + } + } + if (Teuchos::nonnull(dxdp)) { + dxdp = app_->getAdaptSolMgr()->getCurrentDxDp(); + } + + if (adaptData->type==AdaptationType::Topology) { + // This should trigger the nonlinear solver to be rebuilt, which should create new linear + // algebra objects (jac and residual) + auto stepper = integrator.getStepper(); + auto stepper_impl = Teuchos::rcp_dynamic_cast>(stepper); + if (Teuchos::nonnull(stepper_impl)) { + stepper_impl->setDefaultSolver(); + } else { + throw std::runtime_error("Cannot handle explicit case yet, sorry.\n"); + } + stepper->initialize(); + } + } + observeSolutionImpl (state_nc->getX(),state_nc->getXDot(),state_nc->getXDotDot(),dxdp,time); +} + +} // namespace Albany diff --git a/src/Albany_PiroTempusObserver.hpp b/src/Albany_PiroTempusObserver.hpp new file mode 100644 index 0000000000..6a46f4c96f --- /dev/null +++ b/src/Albany_PiroTempusObserver.hpp @@ -0,0 +1,31 @@ +//*****************************************************************// +// Albany 3.0: Copyright 2016 Sandia Corporation // +// This Software is released under the BSD license detailed // +// in the file "license.txt" in the top-level Albany directory // +//*****************************************************************// + +#ifndef ALBANY_PIRO_TEMPUS_OBSERVER_HPP +#define ALBANY_PIRO_TEMPUS_OBSERVER_HPP + +#include "Albany_PiroObserver.hpp" +#include "Tempus_IntegratorObserverBasic.hpp" + +namespace Albany { + +class PiroTempusObserver : public PiroObserver, + public Tempus::IntegratorObserverBasic +{ +public: + PiroTempusObserver(const Teuchos::RCP& app, + const Teuchos::RCP& model); + + // Observe the end of each time step in the time loop + void observeEndTimeStep(const Tempus::Integrator& integrator) override; +protected: + + Teuchos::RCP app_; +}; + +} // namespace Albany + +#endif // ALBANY_PIRO_TEMPUS_OBSERVER_HPP diff --git a/src/Albany_SolverFactory.cpp b/src/Albany_SolverFactory.cpp index 28dc854e20..7429473d4d 100644 --- a/src/Albany_SolverFactory.cpp +++ b/src/Albany_SolverFactory.cpp @@ -6,6 +6,7 @@ #include "Albany_SolverFactory.hpp" #include "Albany_PiroObserver.hpp" +#include "Albany_PiroTempusObserver.hpp" #include "Albany_ModelEvaluator.hpp" #include "Albany_Application.hpp" #include "Albany_Utils.hpp" @@ -255,7 +256,13 @@ createSolver (const Teuchos::RCP& model_tmp, const auto app = model_tmp->getAlbanyApp(); - auto observer = Teuchos::rcp(new PiroObserver(app, modelWithSolve)); + Teuchos::RCP observer; + + if (solutionMethod=="Transient") { + observer = Teuchos::rcp(new PiroTempusObserver(app, modelWithSolve)); + } else { + observer = Teuchos::rcp(new PiroObserver(app, modelWithSolve)); + } Piro::SolverFactory piroFactory; return piroFactory.createSolver( diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index f190778ee8..f212ba297a 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -67,6 +67,7 @@ set (SOURCES Albany_NullSpaceUtils.cpp Albany_ObserverImpl.cpp Albany_PiroObserver.cpp + Albany_PiroTempusObserver.cpp Albany_StatelessObserverImpl.cpp Albany_StateManager.cpp PHAL_Utilities.cpp From bef80b36c85ed19964486715761accb3366bfb67 Mon Sep 17 00:00:00 2001 From: Luca Bertagna Date: Wed, 8 Jan 2025 16:46:23 -0700 Subject: [PATCH 08/13] Add dummy problem to check that the adaptivity hooks work * Add dummy adaptation strategy for STK disc and 1d problems * Add new adv-diff test with adaptation --- src/disc/stk/Albany_GenericSTKMeshStruct.cpp | 4 + src/disc/stk/Albany_GenericSTKMeshStruct.hpp | 2 + src/disc/stk/Albany_STKDiscretization.cpp | 115 +++++++++++++++- tests/corePDEs/AdvDiff/CMakeLists.txt | 6 + tests/corePDEs/AdvDiff/input_adapt.yaml | 135 +++++++++++++++++++ 5 files changed, 255 insertions(+), 7 deletions(-) create mode 100644 tests/corePDEs/AdvDiff/input_adapt.yaml diff --git a/src/disc/stk/Albany_GenericSTKMeshStruct.cpp b/src/disc/stk/Albany_GenericSTKMeshStruct.cpp index 26dae5f912..b9d692a8d6 100644 --- a/src/disc/stk/Albany_GenericSTKMeshStruct.cpp +++ b/src/disc/stk/Albany_GenericSTKMeshStruct.cpp @@ -158,6 +158,9 @@ setFieldData (const Teuchos::RCP& comm, writeCoordsToMMFile = params->get("Write Coordinates to MatrixMarket", false); transferSolutionToCoords = params->get("Transfer Solution to Coordinates", false); + + // Store a copy of the state info struct + sis_ = sis; } void GenericSTKMeshStruct::setAllPartsIO() @@ -1375,6 +1378,7 @@ GenericSTKMeshStruct::getValidGenericSTKParameters(std::string listname) const validPL->set("Sensitivity Method", "None", "Type of sensitivities requested"); validPL->set("Response Function Index", 0, "Response function index (for adjoint transient sensitivities)"); validPL->set("Sensitivity Parameter Index", 0, "Parameter sensitivity index (for transient sensitivities)"); + validPL->sublist("Mesh Adaptivity", false, "Sublist containing mesh adaptivity options"); return validPL; } diff --git a/src/disc/stk/Albany_GenericSTKMeshStruct.hpp b/src/disc/stk/Albany_GenericSTKMeshStruct.hpp index 20a4f74c08..338bd84cfd 100644 --- a/src/disc/stk/Albany_GenericSTKMeshStruct.hpp +++ b/src/disc/stk/Albany_GenericSTKMeshStruct.hpp @@ -126,6 +126,8 @@ class GenericSTKMeshStruct : public AbstractSTKMeshStruct std::vector m_nodesets_from_sidesets; int num_params; + + Teuchos::RCP sis_; }; } // namespace Albany diff --git a/src/disc/stk/Albany_STKDiscretization.cpp b/src/disc/stk/Albany_STKDiscretization.cpp index 42741461f0..8995bd4fd8 100644 --- a/src/disc/stk/Albany_STKDiscretization.cpp +++ b/src/disc/stk/Albany_STKDiscretization.cpp @@ -14,6 +14,7 @@ #include "Albany_Utils.hpp" #include "Albany_GlobalLocalIndexer.hpp" #include "STKConnManager.hpp" +#include "Albany_TmplSTKMeshStruct.hpp" #include #include @@ -1015,6 +1016,7 @@ STKDiscretization::getSolutionDxDp( solutionFieldContainer->fillSolnSensitivity(result, getDOFManager(), overlapped); } + /*****************************************************************/ /*** Private functions follow. These are just used in above code */ /*****************************************************************/ @@ -2444,18 +2446,117 @@ create_dof_mgr (const std::string& part_name, Teuchos::RCP STKDiscretization:: -checkForAdaptation (const Teuchos::RCP& /* solution */, - const Teuchos::RCP& /* solution_dot */, - const Teuchos::RCP& /* solution_dotdot */, - const Teuchos::RCP& /* dxdp */) const +checkForAdaptation (const Teuchos::RCP& solution, + const Teuchos::RCP& solution_dot, + const Teuchos::RCP& solution_dotdot, + const Teuchos::RCP& dxdp) const { - return Teuchos::rcp(new AdaptationData()); + auto adapt_data = Teuchos::rcp(new AdaptationData()); + + // Only do adaptation for simple 1d problems + auto mesh1d = Teuchos::rcp_dynamic_cast>(stkMeshStruct); + if (mesh1d.is_null()) { + std::cout << "NOT a STK1D mesh...\n"; + return adapt_data; + } + auto& adapt_params = discParams->sublist("Mesh Adaptivity"); + auto adapt_type = adapt_params.get("Type","None"); + if (adapt_type=="None") { + std::cout << "NO adaptation...\n"; + return adapt_data; + } + TEUCHOS_TEST_FOR_EXCEPTION (adapt_type!="Uniform-Hessian-Based", std::runtime_error, + "Error! Adaptation type '" << adapt_type << "' not supported.\n" + " - valid choices: None, Uniform-Hessian-Based\n"); + + double tol = adapt_params.get("Max Hessian"); + auto data = getLocalData(solution); + // Simple check: refine if a proxy of the hessian of x is larger than a tolerance + // TODO: replace with + // 1. if |C_i| > threshold, mark for refinement the whole mesh + // 2. Interpolate solution (and all elem/node fields if possible, but not necessary for adv-diff example) + int num_nodes = data.size(); + std::cout << "adapt mesh, curr sol:"; + for (int i=0; ix = solution; + adapt_data->x_dot = solution_dot; + adapt_data->x_dotdot = solution_dotdot; + adapt_data->dxdp = dxdp; + for (int i=1; itol) { + printf("hess=%f, hprev=%f, hnext=%f,up=%f, uc=%f, un=%f\n",hess,h_prev,h_next,data[i-1],data[i],data[i+1]); + adapt_data->type = AdaptationType::Topology; + break; + } + } + + return adapt_data; } void STKDiscretization:: -adapt (const Teuchos::RCP& /* adaptData */) +adapt (const Teuchos::RCP& adaptData) { - throw NotYetImplemented ("STKDiscretization::adapt"); + // Not sure if we allow calling adapt in general, but just in case + if (adaptData->type==AdaptationType::None) { + return; + } + + TEUCHOS_TEST_FOR_EXCEPTION (adaptData->type!=AdaptationType::Topology, std::runtime_error, + "Error! Adaptation type not supported. Only 'None' and 'Topology' are currently supported.\n"); + + // Solution oscillates. We need to half dx + std::cout << "adapting mesh:\n" + << " creating new 1d mesh..." << std::flush;; + auto mesh1d = Teuchos::rcp_dynamic_cast>(stkMeshStruct); + int num_params = mesh1d->getNumParams(); + int ne_x = discParams->get("1D Elements"); + auto& adapt_params = discParams->sublist("Mesh Adaptivity"); + int factor = adapt_params.get("Refining Factor",2); + discParams->set("1D Elements",factor*ne_x); + stkMeshStruct = Teuchos::rcp(new TmplSTKMeshStruct<1>(discParams,comm,num_params)); + stkMeshStruct->setFieldData(comm,mesh1d->sis_); + this->setFieldData(mesh1d->sis_); + stkMeshStruct->setBulkData(comm); + std::cout << "done!\n" << std::flush;; + + std::cout << " old ne_x: " << ne_x << "\n"; + std::cout << " new ne_x: " << factor*ne_x << "\n"; + + updateMesh(); + + int num_time_deriv = discParams->get("Number Of Time Derivatives"); + auto x_mv_new = Thyra::createMembers(getVectorSpace(),num_time_deriv); + + for (int ideriv=0; iderivcol(ideriv)); + auto x = ideriv==0 ? adaptData->x : (ideriv==1 ? adaptData->x_dot : adaptData->x_dotdot); + auto data_old = getLocalData(x); + int num_nodes_new = data_new.size(); + + std::cout << " new solution space size: " << num_nodes_new << "\n"; + for (int inode=0; inode(rem) / factor; + data_new[inode] = data_old[coarse]*(1-alpha) + data_old[coarse+1]*alpha; + } + } + } + + writeSolutionMVToMeshDatabase(*x_mv_new, Teuchos::null, 0, false); } } // namespace Albany diff --git a/tests/corePDEs/AdvDiff/CMakeLists.txt b/tests/corePDEs/AdvDiff/CMakeLists.txt index 9fff09f742..fd4540c810 100644 --- a/tests/corePDEs/AdvDiff/CMakeLists.txt +++ b/tests/corePDEs/AdvDiff/CMakeLists.txt @@ -10,4 +10,10 @@ IF(NOT ALBANY_PARALLEL_ONLY AND ALBANY_IFPACK2) ${CMAKE_CURRENT_BINARY_DIR}/input.yaml COPYONLY) add_test(${testName} ${SerialAlbany.exe} input.yaml) set_tests_properties(${testName} PROPERTIES LABELS "Demo;Forward;Serial") + + # Create adaptivity test + configure_file(${CMAKE_CURRENT_SOURCE_DIR}/input_adapt.yaml + ${CMAKE_CURRENT_BINARY_DIR}/input_adapt.yaml COPYONLY) + add_test(${testName}_adapt ${SerialAlbany.exe} input_adapt.yaml) + set_tests_properties(${testName} PROPERTIES LABELS "Demo;Forward;Serial;Adapt") ENDIF () diff --git a/tests/corePDEs/AdvDiff/input_adapt.yaml b/tests/corePDEs/AdvDiff/input_adapt.yaml new file mode 100644 index 0000000000..081895b7e8 --- /dev/null +++ b/tests/corePDEs/AdvDiff/input_adapt.yaml @@ -0,0 +1,135 @@ +%YAML 1.1 +--- +ANONYMOUS: + Debug Output: + Write Solution to MatrixMarket: 0 + Report Timers: false + Problem: + Name: AdvDiff 1D + Solution Method: Transient + Number of PDE Equations: 1 + Dirichlet BCs: + DBC on NS NodeSet0 for DOF U0: 0.0 + DBC on NS NodeSet1 for DOF U0: 1.0 + Initial Condition: + Function: Constant + Function Data: [0.0] + Options: + Advection Coefficient: [10.0] + Viscosity mu: 0.5 + Parameters: + Number Of Parameters: 0 + Response Functions: + Number Of Responses: 3 + Response 2: + Type: Scalar Response + Name: Solution Average + Response 1: + Equation: 1 + Type: Scalar Response + Name: Solution Max Value + Response 0: + Equation: 0 + Type: Scalar Response + Name: Solution Max Value + Discretization: + 1D Elements: 10 + 1D Scale: 1.0 + Method: STK1D + Exodus Output File Name: adv_diff1D_adapt.exo + Mesh Adaptivity: + Type: Uniform-Hessian-Based + Max Hessian: 1.0 + Refining Factor: 2 + Piro: + Tempus: + Integrator Name: Tempus Integrator + Tempus Integrator: + Integrator Type: Integrator Basic + Stepper Name: Tempus Stepper + Solution History: + Storage Type: Unlimited + Storage Limit: 20 + Time Step Control: + Initial Time: 0.0 + Initial Time Index: 0 + Initial Time Step: 0.025 + Final Time: 1.0 + Maximum Absolute Error: 1.0e-08 + Maximum Relative Error: 1.0e-08 + Maximum Number of Stepper Failures: 10 + Maximum Number of Consecutive Stepper Failures: 5 + Tempus Stepper: + Stepper Type: Backward Euler + Solver Name: Demo Solver + Demo Solver: + NOX: + Direction: + Method: Newton + Newton: + Forcing Term Method: Constant + Rescue Bad Newton Solve: true + Linear Solver: + Tolerance: 1.0e-05 + Line Search: + Full Step: + Full Step: 1.0 + Method: Full Step + Nonlinear Solver: Line Search Based + Printing: + Output Precision: 3 + Output Processor: 0 + Output Information: + Error: true + Warning: true + Outer Iteration: false + Parameters: true + Details: false + Linear Solver Details: true + Stepper Iteration: true + Stepper Details: true + Stepper Parameters: true + Solver Options: + Status Test Check Type: Minimal + Status Tests: + Test Type: Combo + Combo Type: OR + Number of Tests: 2 + Test 0: + Test Type: NormF + Tolerance: 1.0e-08 + Test 1: + Test Type: MaxIters + Maximum Iterations: 10 + Stratimikos: + Linear Solver Type: Belos + Linear Solver Types: + Belos: + Solver Type: Block GMRES + Solver Types: + Block GMRES: + Convergence Tolerance: 1.0e-05 + Output Frequency: 10 + Output Style: 1 + Verbosity: 0 + Maximum Iterations: 10 + Block Size: 1 + Num Blocks: 100 + Flexible Gmres: false + Preconditioner Type: Ifpack2 + Preconditioner Types: + Ifpack2: + Prec Type: ILUT + Overlap: 1 + Ifpack2 Settings: + 'fact: ilut level-of-fill': 1.0e+00 + Regression For Response 0: + Test Value: 2.435402154398e-01 + Relative Tolerance: 1.0e-04 + Regression For Response 1: + Test Value: 2.915474759851e-01 + Relative Tolerance: 1.0e-04 + Regression For Response 2: + Test Value: 5.749999999897e-02 + Relative Tolerance: 1.0e-04 +... From 6b711ef7c78d2aa41f9853a374053842c88c1b5f Mon Sep 17 00:00:00 2001 From: Luca Bertagna Date: Wed, 8 Jan 2025 18:34:13 -0700 Subject: [PATCH 09/13] Some fixes related to const-ness --- src/Albany_PiroTempusObserver.cpp | 6 ++++-- src/SolutionManager.hpp | 6 ++++-- 2 files changed, 8 insertions(+), 4 deletions(-) diff --git a/src/Albany_PiroTempusObserver.cpp b/src/Albany_PiroTempusObserver.cpp index e8b6767b00..b059ba3f71 100644 --- a/src/Albany_PiroTempusObserver.cpp +++ b/src/Albany_PiroTempusObserver.cpp @@ -50,7 +50,9 @@ observeEndTimeStep(const Tempus::Integrator& integrator) auto time = state_nc->getTime(); auto disc = app_->getDiscretization(); - auto adaptData = disc->checkForAdaptation(state_nc->getX(),state_nc->getXDot(),state_nc->getXDotDot(),dxdp); + + auto state = state_nc.getConst(); + auto adaptData = disc->checkForAdaptation(state->getX(),state->getXDot(),state->getXDotDot(),dxdp); // Before observing the solution, check if we need to adapt if (adaptData->type!=AdaptationType::None) { @@ -87,7 +89,7 @@ observeEndTimeStep(const Tempus::Integrator& integrator) stepper->initialize(); } } - observeSolutionImpl (state_nc->getX(),state_nc->getXDot(),state_nc->getXDotDot(),dxdp,time); + observeSolutionImpl (state->getX(),state->getXDot(),state->getXDotDot(),dxdp,time); } } // namespace Albany diff --git a/src/SolutionManager.hpp b/src/SolutionManager.hpp index 1927382dcb..b4862dc9a1 100644 --- a/src/SolutionManager.hpp +++ b/src/SolutionManager.hpp @@ -42,8 +42,10 @@ class SolutionManager { Teuchos::RCP get_cas_manager() const { return cas_manager; } - Teuchos::RCP getCurrentSolution() { return current_soln; } - Teuchos::RCP getCurrentDxDp() { return current_dxdp; } + Teuchos::RCP getCurrentSolution() { return current_soln; } + Teuchos::RCP getCurrentDxDp() { return current_dxdp; } + Teuchos::RCP getCurrentSolution() const { return current_soln; } + Teuchos::RCP getCurrentDxDp() const { return current_dxdp; } void scatterX (const Thyra_Vector& x, const Teuchos::Ptr x_dot, From e2c9dc9579a83696ad8e755ed6e4e06434e8b798 Mon Sep 17 00:00:00 2001 From: Luca Bertagna Date: Wed, 8 Jan 2025 18:35:07 -0700 Subject: [PATCH 10/13] Better handle the rebuilding of the tempus stepper nonlinear solver during adaptation --- src/Albany_PiroTempusObserver.cpp | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/src/Albany_PiroTempusObserver.cpp b/src/Albany_PiroTempusObserver.cpp index b059ba3f71..0e8c48ed80 100644 --- a/src/Albany_PiroTempusObserver.cpp +++ b/src/Albany_PiroTempusObserver.cpp @@ -80,12 +80,7 @@ observeEndTimeStep(const Tempus::Integrator& integrator) // This should trigger the nonlinear solver to be rebuilt, which should create new linear // algebra objects (jac and residual) auto stepper = integrator.getStepper(); - auto stepper_impl = Teuchos::rcp_dynamic_cast>(stepper); - if (Teuchos::nonnull(stepper_impl)) { - stepper_impl->setDefaultSolver(); - } else { - throw std::runtime_error("Cannot handle explicit case yet, sorry.\n"); - } + stepper->setModel(model_); stepper->initialize(); } } From 2eb4abff2bcb4d7bc4daa8ad268d5c0f09af020f Mon Sep 17 00:00:00 2001 From: Luca Bertagna Date: Fri, 10 Jan 2025 13:44:18 -0700 Subject: [PATCH 11/13] Fix adaptation unit test --- src/disc/stk/Albany_STKDiscretization.cpp | 23 ++++++----------------- tests/corePDEs/AdvDiff/CMakeLists.txt | 2 +- tests/corePDEs/AdvDiff/input_adapt.yaml | 22 ++++++++++------------ 3 files changed, 17 insertions(+), 30 deletions(-) diff --git a/src/disc/stk/Albany_STKDiscretization.cpp b/src/disc/stk/Albany_STKDiscretization.cpp index 8995bd4fd8..037e9f494a 100644 --- a/src/disc/stk/Albany_STKDiscretization.cpp +++ b/src/disc/stk/Albany_STKDiscretization.cpp @@ -2462,12 +2462,11 @@ checkForAdaptation (const Teuchos::RCP& solution, auto& adapt_params = discParams->sublist("Mesh Adaptivity"); auto adapt_type = adapt_params.get("Type","None"); if (adapt_type=="None") { - std::cout << "NO adaptation...\n"; return adapt_data; } - TEUCHOS_TEST_FOR_EXCEPTION (adapt_type!="Uniform-Hessian-Based", std::runtime_error, + TEUCHOS_TEST_FOR_EXCEPTION (adapt_type!="Minimally-Oscillatory", std::runtime_error, "Error! Adaptation type '" << adapt_type << "' not supported.\n" - " - valid choices: None, Uniform-Hessian-Based\n"); + " - valid choices: None, Minimally-Oscillatory\n"); double tol = adapt_params.get("Max Hessian"); auto data = getLocalData(solution); @@ -2476,11 +2475,6 @@ checkForAdaptation (const Teuchos::RCP& solution, // 1. if |C_i| > threshold, mark for refinement the whole mesh // 2. Interpolate solution (and all elem/node fields if possible, but not necessary for adv-diff example) int num_nodes = data.size(); - std::cout << "adapt mesh, curr sol:"; - for (int i=0; ix = solution; @@ -2491,8 +2485,9 @@ checkForAdaptation (const Teuchos::RCP& solution, auto h_prev = coordinates[i] - coordinates[i-1]; auto h_next = coordinates[i+1] - coordinates[i]; auto hess = (data[i-1] - 2*data[i] + data[i+1]) / (h_prev*h_next); - if (std::fabs(hess)>tol) { - printf("hess=%f, hprev=%f, hnext=%f,up=%f, uc=%f, un=%f\n",hess,h_prev,h_next,data[i-1],data[i],data[i+1]); + auto grad_prev = (data[i]-data[i-1]) / h_prev; + auto grad_next = (data[i+1]-data[i]) / h_next; + if (std::fabs(hess)>tol and grad_prev*grad_next<0) { adapt_data->type = AdaptationType::Topology; break; } @@ -2513,22 +2508,17 @@ adapt (const Teuchos::RCP& adaptData) "Error! Adaptation type not supported. Only 'None' and 'Topology' are currently supported.\n"); // Solution oscillates. We need to half dx - std::cout << "adapting mesh:\n" - << " creating new 1d mesh..." << std::flush;; auto mesh1d = Teuchos::rcp_dynamic_cast>(stkMeshStruct); int num_params = mesh1d->getNumParams(); int ne_x = discParams->get("1D Elements"); auto& adapt_params = discParams->sublist("Mesh Adaptivity"); + discParams->set("Workset Size", stkMeshStruct->meshSpecs()[0]->worksetSize); int factor = adapt_params.get("Refining Factor",2); discParams->set("1D Elements",factor*ne_x); stkMeshStruct = Teuchos::rcp(new TmplSTKMeshStruct<1>(discParams,comm,num_params)); stkMeshStruct->setFieldData(comm,mesh1d->sis_); this->setFieldData(mesh1d->sis_); stkMeshStruct->setBulkData(comm); - std::cout << "done!\n" << std::flush;; - - std::cout << " old ne_x: " << ne_x << "\n"; - std::cout << " new ne_x: " << factor*ne_x << "\n"; updateMesh(); @@ -2541,7 +2531,6 @@ adapt (const Teuchos::RCP& adaptData) auto data_old = getLocalData(x); int num_nodes_new = data_new.size(); - std::cout << " new solution space size: " << num_nodes_new << "\n"; for (int inode=0; inode Date: Wed, 12 Mar 2025 18:00:24 -0600 Subject: [PATCH 12/13] Remove stranded debug line --- src/disc/stk/Albany_STKDiscretization.cpp | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/src/disc/stk/Albany_STKDiscretization.cpp b/src/disc/stk/Albany_STKDiscretization.cpp index 037e9f494a..5ecdaaf596 100644 --- a/src/disc/stk/Albany_STKDiscretization.cpp +++ b/src/disc/stk/Albany_STKDiscretization.cpp @@ -2453,12 +2453,6 @@ checkForAdaptation (const Teuchos::RCP& solution, { auto adapt_data = Teuchos::rcp(new AdaptationData()); - // Only do adaptation for simple 1d problems - auto mesh1d = Teuchos::rcp_dynamic_cast>(stkMeshStruct); - if (mesh1d.is_null()) { - std::cout << "NOT a STK1D mesh...\n"; - return adapt_data; - } auto& adapt_params = discParams->sublist("Mesh Adaptivity"); auto adapt_type = adapt_params.get("Type","None"); if (adapt_type=="None") { @@ -2468,6 +2462,11 @@ checkForAdaptation (const Teuchos::RCP& solution, "Error! Adaptation type '" << adapt_type << "' not supported.\n" " - valid choices: None, Minimally-Oscillatory\n"); + // Only do adaptation for simple 1d problems + auto mesh1d = Teuchos::rcp_dynamic_cast>(stkMeshStruct); + TEUCHOS_TEST_FOR_EXCEPTION (mesh1d.is_null(), std::runtime_error, + "Error! Adaptation for STK is only supported for a simple 1D problem, with STK1D discretization.\n"); + double tol = adapt_params.get("Max Hessian"); auto data = getLocalData(solution); // Simple check: refine if a proxy of the hessian of x is larger than a tolerance From 0895ff5cacc095ef5b88d20281bc298465c80203 Mon Sep 17 00:00:00 2001 From: Luca Bertagna Date: Wed, 12 Mar 2025 18:00:58 -0600 Subject: [PATCH 13/13] Fix handling of Tempus state that stores sensitivities --- src/Albany_PiroTempusObserver.cpp | 45 ++++++++++++++++++++++++------- 1 file changed, 35 insertions(+), 10 deletions(-) diff --git a/src/Albany_PiroTempusObserver.cpp b/src/Albany_PiroTempusObserver.cpp index 0e8c48ed80..2166780bcc 100644 --- a/src/Albany_PiroTempusObserver.cpp +++ b/src/Albany_PiroTempusObserver.cpp @@ -6,7 +6,8 @@ #include "Albany_PiroTempusObserver.hpp" -#include +#include + #include #include @@ -43,16 +44,36 @@ observeEndTimeStep(const Tempus::Integrator& integrator) // since we don't want to write the solution to file, or to observe responses Teuchos::RCP dxdp; auto integrator_ptr = Teuchos::rcpFromRef(integrator); - auto fwd_sens_integrator = Teuchos::rcp_dynamic_cast>(integrator_ptr); - if (Teuchos::nonnull(fwd_sens_integrator)) { - dxdp = fwd_sens_integrator->getDxDp(); - } auto time = state_nc->getTime(); auto disc = app_->getDiscretization(); auto state = state_nc.getConst(); - auto adaptData = disc->checkForAdaptation(state->getX(),state->getXDot(),state->getXDotDot(),dxdp); + auto x = state->getX(); + auto xdot = state->getXDot(); + auto xdotdot = state->getXDotDot(); + + // If piro created a sensitivity tempus integrator, x/xdot/xdotdot will + // be in fact product vectors, so we need to extract the pieces + using DMVPV = Thyra::DefaultMultiVectorProductVector; + auto px = Teuchos::rcp_dynamic_cast(x); + auto pxdot = Teuchos::rcp_dynamic_cast(xdot); + auto pxdotdot = Teuchos::rcp_dynamic_cast(xdotdot); + if (Teuchos::nonnull(px)) { + x = px->getMultiVector()->col(0); + if (Teuchos::nonnull(pxdot)) { + xdot = pxdot->getMultiVector()->col(0); + if (Teuchos::nonnull(pxdotdot)) { + xdotdot = pxdotdot->getMultiVector()->col(0); + } + } + + const int num_param = px->getMultiVector()->domain()->dim() - 1; + const Teuchos::Range1D rng(1, num_param); + dxdp = px->getMultiVector()->subView(rng); + } + + auto adaptData = disc->checkForAdaptation(x,xdot,xdotdot,dxdp); // Before observing the solution, check if we need to adapt if (adaptData->type!=AdaptationType::None) { @@ -65,11 +86,15 @@ observeEndTimeStep(const Tempus::Integrator& integrator) auto sol = app_->getAdaptSolMgr()->getCurrentSolution(); // Reset vectors now that they have been adapted - state_nc->setX(sol->col(0)); + // TODO: we may need to revise these lines in case the state stores product vectors. + x = sol->col(0); + state_nc->setX(x); if (num_time_derivs>0) { - state_nc->setXDot(sol->col(1)); + xdot = sol->col(1); + state_nc->setXDot(xdot); if (num_time_derivs>1) { - state_nc->setXDotDot(sol->col(2)); + xdotdot = sol->col(2); + state_nc->setXDotDot(xdotdot); } } if (Teuchos::nonnull(dxdp)) { @@ -84,7 +109,7 @@ observeEndTimeStep(const Tempus::Integrator& integrator) stepper->initialize(); } } - observeSolutionImpl (state->getX(),state->getXDot(),state->getXDotDot(),dxdp,time); + observeSolutionImpl (x,xdot,xdotdot,dxdp,time); } } // namespace Albany