From c1f6a42c42933dbbf41687fafd9faaf0a550ad7c Mon Sep 17 00:00:00 2001 From: Roger Pawlowski Date: Fri, 20 Sep 2024 15:34:28 -0600 Subject: [PATCH] Panzer: added ROL transient optimization problem Signed-off-by: Roger Pawlowski --- .../example/main_driver/CMakeLists.txt | 6 +- .../energy-transient-tempus-opt-blocked.xml | 204 +++++++- .../example/main_driver/main_driver_opt.cpp | 331 ++----------- .../example/main_driver/rol_params.xml | 109 +---- .../user_app_TempusObserver_WriteToExodus.hpp | 2 + .../user_app_TransientObjective.cpp | 5 + .../user_app_TransientObjective.hpp | 91 ++++ .../user_app_TransientObjective_impl.hpp | 455 ++++++++++++++++++ .../main_driver/user_app_Utilities.cpp | 309 ++++++++++++ .../main_driver/user_app_Utilities.hpp | 145 ++++++ .../src/Panzer_ModelEvaluator_impl.hpp | 10 +- .../src/evaluators/Panzer_Parameter_impl.hpp | 3 +- .../user_app_EquationSetFactory.hpp | 5 + 13 files changed, 1265 insertions(+), 410 deletions(-) create mode 100644 packages/panzer/adapters-stk/example/main_driver/user_app_TransientObjective.cpp create mode 100644 packages/panzer/adapters-stk/example/main_driver/user_app_TransientObjective.hpp create mode 100644 packages/panzer/adapters-stk/example/main_driver/user_app_TransientObjective_impl.hpp create mode 100644 packages/panzer/adapters-stk/example/main_driver/user_app_Utilities.cpp create mode 100644 packages/panzer/adapters-stk/example/main_driver/user_app_Utilities.hpp diff --git a/packages/panzer/adapters-stk/example/main_driver/CMakeLists.txt b/packages/panzer/adapters-stk/example/main_driver/CMakeLists.txt index 5dff70456774..9de8071715b0 100644 --- a/packages/panzer/adapters-stk/example/main_driver/CMakeLists.txt +++ b/packages/panzer/adapters-stk/example/main_driver/CMakeLists.txt @@ -117,11 +117,11 @@ IF(${PACKAGE_NAME}_ENABLE_Teko) PASS_REGULAR_EXPRESSION "panzer::MainDriver run completed." ) - # Optimization with ROL - IF(${PACKAGE_NAME}_ENABLE_ROL) + # Optimization with ROL (currently doesn't support GPU devices) + IF(${PACKAGE_NAME}_ENABLE_ROL AND NOT Kokkos_ENABLE_CUDA AND NOT Kokkos_ENABLE_HIP) TRIBITS_ADD_EXECUTABLE( main_driver_opt - SOURCES main_driver_opt.cpp ${main_driver_SOURCES} user_app_ROLTempusReducedObjective.hpp + SOURCES main_driver_opt.cpp ${main_driver_SOURCES} user_app_ROLTempusReducedObjective.hpp user_app_TransientObjective.hpp user_app_TransientObjective_impl.hpp user_app_TransientObjective.cpp user_app_Utilities.hpp user_app_Utilities.cpp ) TRIBITS_ADD_ADVANCED_TEST( main_driver_energy-transient-tempus-opt-blocked diff --git a/packages/panzer/adapters-stk/example/main_driver/energy-transient-tempus-opt-blocked.xml b/packages/panzer/adapters-stk/example/main_driver/energy-transient-tempus-opt-blocked.xml index fc06b2ef718d..5e377070a278 100644 --- a/packages/panzer/adapters-stk/example/main_driver/energy-transient-tempus-opt-blocked.xml +++ b/packages/panzer/adapters-stk/example/main_driver/energy-transient-tempus-opt-blocked.xml @@ -288,7 +288,7 @@ - + @@ -394,4 +394,206 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/packages/panzer/adapters-stk/example/main_driver/main_driver_opt.cpp b/packages/panzer/adapters-stk/example/main_driver/main_driver_opt.cpp index af09ef0223cf..c0fae5d462e6 100644 --- a/packages/panzer/adapters-stk/example/main_driver/main_driver_opt.cpp +++ b/packages/panzer/adapters-stk/example/main_driver/main_driver_opt.cpp @@ -19,7 +19,6 @@ #include "Teuchos_XMLParameterListHelpers.hpp" #include "Teuchos_YamlParameterListHelpers.hpp" #include "Teuchos_FancyOStream.hpp" -#include "Teuchos_oblackholestream.hpp" #include "Teuchos_Assert.hpp" #include "Teuchos_as.hpp" @@ -49,6 +48,7 @@ #include "user_app_ResponseEvaluatorFactory_HOFlux.hpp" #include "Panzer_ROLTempusReducedObjective.hpp" +#include "user_app_TransientObjective.hpp" #include "Panzer_ResponseEvaluatorFactory_Probe.hpp" @@ -58,38 +58,16 @@ #include +#include "user_app_Utilities.hpp" + #include #include #include -// ************************************************* -// applicaiton functions defined in this file -// ************************************************* -namespace user_app { - void addResponsesToModelEvaluatorFactory(const Teuchos::ParameterList& response_sublist, - panzer_stk::ModelEvaluatorFactory& me_factory); - - void computeAndPrintResponses(const Teuchos::RCP>& me, - const auto& x, - const auto& x_dot, - const auto& t, - const auto& global_data, - std::ostream& out); - - std::tuple>, - Teuchos::RCP, - Teuchos::RCP>> - buildTempusIntegrator(const Teuchos::RCP& input_params, - const Teuchos::RCP>& comm, - const bool overrideNoxOutput); -} - -// ************************************************* // ************************************************* // Main Driver for a transient optimization problem using ROL and // Tempus // ************************************************* -// ************************************************* int main(int argc, char *argv[]) { @@ -98,8 +76,8 @@ int main(int argc, char *argv[]) int status = 0; - std::ostream blackhole{nullptr}; - Teuchos::GlobalMPISession mpiSession(&argc, &argv, &blackhole); + std::ostream os_dev_null{nullptr}; + Teuchos::GlobalMPISession mpiSession(&argc, &argv, &os_dev_null); Kokkos::initialize(argc,argv); Teuchos::RCP out = Teuchos::rcp(new Teuchos::FancyOStream(Teuchos::rcp(&std::cout,false))); @@ -155,44 +133,43 @@ int main(int argc, char *argv[]) if (printInputPL) *out << *input_params << std::endl; + // Pull off objective and rol lists so the panzer validation does not fail + RCP objective_params = parameterList(input_params->sublist("Objective",true)); + RCP rol_params = parameterList(input_params->sublist("ROL", true)); + input_params->remove("Objective"); + input_params->remove("ROL"); - RCP params_copy = parameterList(*input_params); - auto [integrator,global_data,physics] = user_app::buildTempusIntegrator(params_copy,comm,overrideNoxOutput); - - bool integratorStatus = integrator->advanceTime(); - TEUCHOS_ASSERT(integratorStatus); - - auto x = integrator->getSolutionHistory()->getCurrentState()->getX(); - auto x_dot = integrator->getSolutionHistory()->getCurrentState()->getXDot(); - auto t = integrator->getSolutionHistory()->getCurrentState()->getTime(); - - user_app::computeAndPrintResponses(physics,x,x_dot,t,global_data,*out); - { - Teuchos::RCP rol_params = Teuchos::rcp(new Teuchos::ParameterList("ROL Parameters")); - { - const std::string rol_input_file = "rol_params.xml"; - Teuchos::updateParametersFromXmlFileAndBroadcast(rol_input_file, rol_params.ptr(), *comm); - - RCP pl = sublist(rol_params, "Tempus", true); - RCP objective_params = sublist(rol_params, "Objective", true); - ROL::Ptr > objective = - ROL::makePtr >(physics, pl, objective_params); - - // Create target -- do forward integration with perturbed parameter values - RCP > zs = objective->create_design_vector(); - zs->setScalar(1.5); - RCP > r = objective->create_response_vector(); - objective->run_tempus(*r, *zs); - objective->set_target(r); - - } - + RCP tempus_params = parameterList(input_params->sublist("Solution Control",true).sublist("Tempus",true)); + auto objective = ROL::makePtr>(input_params,comm,objective_params,out); + + // Create target -- do forward integration with perturbed parameter values + RCP> p = objective->create_design_vector(); + p->setScalar(2.0); + RCP> r = objective->create_response_vector(); + objective->run_tempus(*r, *p); + objective->set_target(r); + + p->setScalar(1.5); + ROL::OptimizationProblem problem(objective, p); + ROL::OptimizationSolver solver(problem, *rol_params); ROL::Ptr outStream = ROL::makePtrFromRef(std::cout); + if (comm->getRank() == 0) + solver.solve(*outStream); + else + solver.solve(os_dev_null); - + { + const ROL::ThyraVector& thyra_p = dyn_cast >(*p); + ROL::ThyraVector& thyra_r = dyn_cast >(*r); + *out << "Final Values: p = " << Thyra::get_ele(*(thyra_p.getVector()),0) + << ", g = " << Thyra::get_ele(*(thyra_r.getVector()),0) + << std::endl; + } } + // user_app::computeAndPrintResponses(physics,x,x_dot,t,global_data,*out); + stackedTimer->stopBaseTimer(); if (printTimers) { Teuchos::StackedTimer::OutputOptions options; @@ -218,239 +195,3 @@ int main(int argc, char *argv[]) return status; } - -// ************************************************* -// ************************************************* - -void user_app::addResponsesToModelEvaluatorFactory(const Teuchos::ParameterList& responses, - panzer_stk::ModelEvaluatorFactory& me_factory) -{ - for(Teuchos::ParameterList::ConstIterator itr=responses.begin();itr!=responses.end();++itr) { - const std::string name = responses.name(itr); - TEUCHOS_ASSERT(responses.entry(itr).isList()); - Teuchos::ParameterList & lst = Teuchos::getValue(responses.entry(itr)); - - if (lst.get("Type") == "Functional") { - - panzer::FunctionalResponse_Builder builder; - builder.comm = MPI_COMM_WORLD; - builder.cubatureDegree = 2; - builder.requiresCellIntegral = true; - builder.quadPointField = lst.get("Field Name"); - - std::vector eblocks; - panzer::StringTokenizer(eblocks,lst.get("Element Blocks"),",",true); - - std::vector wkst_descs; - for(std::size_t i=0;i("Type") == "Point Value") { - panzer::ProbeResponse_Builder builder; - builder.comm = MPI_COMM_WORLD; - builder.point = Teuchos::Array{0.5,0.5}; - builder.cubatureDegree = 2; - builder.fieldName = lst.get("Field Name"); - builder.applyDirichletToDerivative = false; - - std::vector eblocks; - panzer::StringTokenizer(eblocks,lst.get("Element Blocks"),",",true); - - std::vector descriptors; - for(std::size_t i=0;i("Type") << "\" is not supported!"); - } - } -} - -void user_app::computeAndPrintResponses(const Teuchos::RCP>& physics, - const auto& x, - const auto& x_dot, - const auto& t, - const auto& global_data, - std::ostream& os) -{ - if(physics->Ng()>0) { - - os << "Ng = " << physics->Ng() << std::endl; - os << "Np = " << physics->Np() << std::endl; - Thyra::ModelEvaluatorBase::InArgs respInArgs = physics->createInArgs(); - Thyra::ModelEvaluatorBase::OutArgs respOutArgs = physics->createOutArgs(); - - TEUCHOS_ASSERT(physics->Ng()==respOutArgs.Ng()); - - respInArgs.set_x(x); - respInArgs.set_x_dot(x_dot); - - for(int g=0;g > response = Thyra::createMember(*physics->get_g_space(g)); - respOutArgs.set_g(g,response); - - for (int p = 0; p < physics->Np(); ++p) { - bool derivative_supported = respOutArgs.supports(Thyra::ModelEvaluatorBase::OUT_ARG_DgDp,g,p).supports(Thyra::ModelEvaluatorBase::DERIV_MV_JACOBIAN_FORM); - os << "DgDp(" << g << "," << p << ") supports = " << derivative_supported << std::endl; - if (derivative_supported) { - auto dgdp = physics->create_DgDp_op(g,p); - respOutArgs.set_DgDp(g,p,Thyra::ModelEvaluatorBase::Derivative(dgdp)); - } - } - } - - physics->evalModel(respInArgs, respOutArgs); - - { - auto parameter_library = global_data->pl; - TEUCHOS_ASSERT(nonnull(parameter_library)); - for (auto i=parameter_library->begin();i != parameter_library->end(); ++i) { - os << "Sacado::ParameterLibrary: " << i->first << std::endl; - } - } - - { - auto nominalValues = physics->getNominalValues(); - for (int i=0; i < respInArgs.Np();++i) { - auto p = nominalValues.get_p(i); - auto p_names = physics->get_p_names(i); - os << "ModelEvaluator::NominalValues Parameter Value: \"" << (*p_names)[0] << "\" = " << Thyra::get_ele(*p,0) << std::endl; - } - } - - for (int i=0;i > response = respOutArgs.get_g(i); - TEUCHOS_ASSERT(response!=Teuchos::null); - os << "Response Value: \"" << physics->get_g_names(i)[0] << "\" = " << Thyra::get_ele(*response,0) << std::endl; - for (int j=0; j < respOutArgs.Np(); ++j) { - // os << " dg(" << i << ")/dp(" << j << ") supports(GRAD_FORM) = " << respOutArgs.supports(Thyra::ModelEvaluatorBase::OUT_ARG_DgDp,i,j).supports(Thyra::ModelEvaluatorBase::DERIV_MV_GRADIENT_FORM) << std::endl; - // os << " dg(" << i << ")/dp(" << j << ") supports(JAC_FORM) = " << respOutArgs.supports(Thyra::ModelEvaluatorBase::OUT_ARG_DgDp,i,j).supports(Thyra::ModelEvaluatorBase::DERIV_MV_JACOBIAN_FORM) << std::endl; - } - } - } -} - -std::tuple>, - Teuchos::RCP, - Teuchos::RCP>> -user_app::buildTempusIntegrator(const Teuchos::RCP& input_params, - const Teuchos::RCP>& comm, - const bool overrideNoxOutput) -{ - using namespace Teuchos; - using namespace Tempus; - - Teuchos::ParameterList solver_factories = input_params->sublist("Solver Factories"); - input_params->remove("Solver Factories"); - - // Add in the application specific equation set factory - Teuchos::RCP eqset_factory = Teuchos::rcp(new user_app::MyFactory); - - // Add in the application specific closure model factory - user_app::MyModelFactory_TemplateBuilder cm_builder; - panzer::ClosureModelFactory_TemplateManager cm_factory; - cm_factory.buildObjects(cm_builder); - - // Add in the application specific bc factory - user_app::BCFactory bc_factory; - - // Create the global data - Teuchos::RCP global_data = panzer::createGlobalData(); - input_params->sublist("User Data").set("Comm", comm); // Required for GlobalStatistics - - Teuchos::RCP > stkIOResponseLibrary - = Teuchos::rcp(new panzer::ResponseLibrary()); - - Teuchos::RCP> physics; - Teuchos::RCP> solver; - Teuchos::RCP> rLibrary; - std::vector> physicsBlocks; - Teuchos::RCP> linObjFactory; - Teuchos::RCP globalIndexer; - Teuchos::RCP mesh; - { - Teuchos::ParameterList responses = input_params->sublist("Responses"); - input_params->remove("Responses"); - - panzer_stk::ModelEvaluatorFactory me_factory; - me_factory.setParameterList(input_params); - auto strat_list = Teuchos::parameterList(); - *strat_list = input_params->sublist("Solution Control",true).sublist("Tempus",true).sublist("My Example Stepper",true).sublist("My Example Solver",true).sublist("NOX",true).sublist("Direction",true).sublist("Newton",true).sublist("Stratimikos Linear Solver",true).sublist("Stratimikos",true); - me_factory.setStratimikosList(strat_list); - me_factory.buildObjects(comm,global_data,eqset_factory,bc_factory,cm_factory); - - user_app::addResponsesToModelEvaluatorFactory(responses,me_factory); - me_factory.buildResponses(cm_factory); - - physicsBlocks = me_factory.getPhysicsBlocks(); - physics = me_factory.getPhysicsModelEvaluator(); - rLibrary = me_factory.getResponseLibrary(); - linObjFactory = me_factory.getLinearObjFactory(); - globalIndexer = me_factory.getGlobalIndexer(); - mesh = me_factory.getMesh(); - } - - // setup outputs to mesh on the stkIOResponseLibrary - { - stkIOResponseLibrary->initialize(*rLibrary); - - Teuchos::ParameterList user_data(input_params->sublist("User Data")); - user_data.set("Workset Size",input_params->sublist("Assembly").get("Workset Size")); - - std::vector eBlocks; - mesh->getElementBlockNames(eBlocks); - - panzer_stk::RespFactorySolnWriter_Builder builder; - builder.mesh = mesh; - stkIOResponseLibrary->addResponse("Main Field Output",eBlocks,builder); - - stkIOResponseLibrary->buildResponseEvaluators(physicsBlocks, - cm_factory, - input_params->sublist("Closure Models"), - user_data); - } - - RCP tempus_pl = parameterList(std::string("Tempus")); - *tempus_pl = input_params->sublist("Solution Control").sublist("Tempus"); - - if (overrideNoxOutput) { - auto nox_utils = Teuchos::rcp(new NOX::Utils(NOX::Utils::Error,comm->getRank(),0,3)); - auto print_nonlinear = Teuchos::rcp(new NOX::ObserverPrint(nox_utils,10)); - tempus_pl->sublist("My Example Stepper",true).sublist("My Example Solver",true).sublist("NOX",true).sublist("Solver Options",true).set>("Observer",print_nonlinear); - } - - // Tempus is overriding NOX parameters with its own defaults. Need - // to fix this. It also does not validate because of arbitrary - // solver name. Should validate but use - // disableRecursiveValidation() to avoid errors with arbitrary - // names. - - const bool doInitialization = false; - auto integrator = createIntegratorBasic(tempus_pl, physics, doInitialization); - - RCP noxList = parameterList("Correct NOX Params"); - *noxList = tempus_pl->sublist("My Example Stepper",true).sublist("My Example Solver",true).sublist("NOX",true); - // noxList->print(std::cout); - integrator->getStepper()->getSolver()->setParameterList(noxList); - integrator->initialize(); - - // Setting observers on tempus breaks the screen output of the - // time steps! It replaces IntegratorObserverBasic which handles - // IO. Is this at least documented? - { - RCP> tempus_observers = rcp(new Tempus::IntegratorObserverComposite); - RCP tof = Teuchos::rcp(new user_app::TempusObserverFactory(stkIOResponseLibrary,rLibrary->getWorksetContainer())); - tempus_observers->addObserver(Teuchos::rcp(new Tempus::IntegratorObserverBasic)); - tempus_observers->addObserver(tof->buildTempusObserver(mesh,globalIndexer,linObjFactory)); - integrator->setObserver(tempus_observers); - } - - RCP> x0 = physics->getNominalValues().get_x()->clone_v(); - integrator->initializeSolutionHistory(0.0, x0); - - return {integrator,global_data,physics}; -} diff --git a/packages/panzer/adapters-stk/example/main_driver/rol_params.xml b/packages/panzer/adapters-stk/example/main_driver/rol_params.xml index 8777e117ca47..27ecf899e374 100644 --- a/packages/panzer/adapters-stk/example/main_driver/rol_params.xml +++ b/packages/panzer/adapters-stk/example/main_driver/rol_params.xml @@ -1,114 +1,9 @@ - - - - - - - - - - - + + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/packages/panzer/adapters-stk/example/main_driver/user_app_TempusObserver_WriteToExodus.hpp b/packages/panzer/adapters-stk/example/main_driver/user_app_TempusObserver_WriteToExodus.hpp index 04cec8bceb14..dc965c8dbf0e 100644 --- a/packages/panzer/adapters-stk/example/main_driver/user_app_TempusObserver_WriteToExodus.hpp +++ b/packages/panzer/adapters-stk/example/main_driver/user_app_TempusObserver_WriteToExodus.hpp @@ -15,6 +15,8 @@ #include "Tempus_IntegratorObserver.hpp" #include "Teuchos_RCP.hpp" #include "Teuchos_Assert.hpp" +#include "Teuchos_TimeMonitor.hpp" +#include "Teuchos_StackedTimer.hpp" #include "Panzer_STK_Interface.hpp" #include "Panzer_GlobalIndexer.hpp" diff --git a/packages/panzer/adapters-stk/example/main_driver/user_app_TransientObjective.cpp b/packages/panzer/adapters-stk/example/main_driver/user_app_TransientObjective.cpp new file mode 100644 index 000000000000..56e1bbf8d563 --- /dev/null +++ b/packages/panzer/adapters-stk/example/main_driver/user_app_TransientObjective.cpp @@ -0,0 +1,5 @@ + +#include "user_app_TransientObjective.hpp" +#include "user_app_TransientObjective_impl.hpp" + +template class ROL::TransientReducedObjective; diff --git a/packages/panzer/adapters-stk/example/main_driver/user_app_TransientObjective.hpp b/packages/panzer/adapters-stk/example/main_driver/user_app_TransientObjective.hpp new file mode 100644 index 000000000000..65b041d76222 --- /dev/null +++ b/packages/panzer/adapters-stk/example/main_driver/user_app_TransientObjective.hpp @@ -0,0 +1,91 @@ +// @HEADER +// @HEADER + +#ifndef USER_APP_TRANSIENT_REDUCED_OBJECTIVE_HPP +#define USER_APP_TRANSIENT_REDUCED_OBJECTIVE_HPP + +#include + +#include "Teuchos_RCP.hpp" +#include "Teuchos_ParameterList.hpp" + +#include "Tempus_IntegratorBasic.hpp" +#include "Tempus_IntegratorForwardSensitivity.hpp" +#include "Tempus_IntegratorAdjointSensitivity.hpp" +#include "Tempus_IntegratorPseudoTransientForwardSensitivity.hpp" +#include "Tempus_IntegratorPseudoTransientAdjointSensitivity.hpp" + +#include "Thyra_ModelEvaluator.hpp" +#include "Thyra_DefaultNominalBoundsOverrideModelEvaluator.hpp" +#include "Thyra_VectorStdOps.hpp" + +#include "ROL_Objective.hpp" +#include "ROL_Vector.hpp" +#include "ROL_ThyraVector.hpp" + +#include "user_app_Utilities.hpp" + +namespace ROL { + +template +class TransientReducedObjective : public virtual ROL::Objective { +public: + + TransientReducedObjective( + const Teuchos::RCP& input_params, + const Teuchos::RCP>& comm, + const Teuchos::RCP& objective_params, + const Teuchos::RCP& os); + + virtual ~TransientReducedObjective() {} + + //! Compute value of objective + Real value( const ROL::Vector &x, Real &tol ); + + //! Compute gradient of objective + void gradient( ROL::Vector &g, const ROL::Vector &x, Real &tol ); + + //! Set response target for computing objective + void set_target(const Teuchos::RCP >& target); + void set_target(const Teuchos::RCP >& target); + + //! Helper function to create optimization vector + Teuchos::RCP > create_design_vector() const; + + //! Helper function to create a response vector + Teuchos::RCP > create_response_vector() const; + + //! Helper function to run tempus, computing responses and derivatives + void run_tempus(ROL::Vector& r, const ROL::Vector& p); + void run_tempus(const Thyra::ModelEvaluatorBase::InArgs& inArgs, + const Thyra::ModelEvaluatorBase::OutArgs& outArgs); + +private: + + // Teuchos::RCP > model_; + Teuchos::RCP tempus_params_; + Teuchos::RCP > target_; + std::string objective_type_; + std::string sensitivity_method_; + int param_index_; + int response_index_; + bool use_fd_gradient_; + + // Objects neeeded to rebuild the time integrator since it seems to + // not be able to reset itself + Teuchos::RCP input_params_; + Teuchos::RCP> comm_; + Teuchos::RCP> model_; + Teuchos::RCP global_data_; + Teuchos::RCP mesh_; + Teuchos::RCP> response_library_; + Teuchos::RCP> stk_io_response_library_; + Teuchos::RCP> lin_obj_factory_; + Teuchos::RCP global_indexer_; + Teuchos::RCP os_; + bool print_debug_; +}; + +} + +#endif diff --git a/packages/panzer/adapters-stk/example/main_driver/user_app_TransientObjective_impl.hpp b/packages/panzer/adapters-stk/example/main_driver/user_app_TransientObjective_impl.hpp new file mode 100644 index 000000000000..41e6053f922a --- /dev/null +++ b/packages/panzer/adapters-stk/example/main_driver/user_app_TransientObjective_impl.hpp @@ -0,0 +1,455 @@ +#include "user_app_Utilities.hpp" +#include "Teuchos_StandardParameterEntryValidators.hpp" + +namespace ROL { + +template +TransientReducedObjective:: +TransientReducedObjective(const Teuchos::RCP& input_params, + const Teuchos::RCP>& comm, + const Teuchos::RCP& objective_params, + const Teuchos::RCP& os) : + objective_type_("Sum of Squares"), + sensitivity_method_("Forward"), + param_index_(0), + response_index_(0), + use_fd_gradient_(false), + input_params_(input_params), + comm_(comm), + os_(os), + print_debug_(true) +{ + TEUCHOS_ASSERT(nonnull(input_params)); + TEUCHOS_ASSERT(nonnull(comm)); + TEUCHOS_ASSERT(nonnull(objective_params)); + + Teuchos::ParameterList valid_params; + valid_params.set("Objective Type", + "Sum of Squares", + "How to compute the objective function", + Teuchos::rcp(new Teuchos::StringValidator({"Sum of Squares","Sum"}))); + valid_params.set("Sensitivity Method", + "Forward", + "Algorithm to compute parameter sensitivities", + Teuchos::rcp(new Teuchos::StringValidator({"Forward", + "Adjoint", + "Pseudotransient Forward", + "Pseudotransient Adjoint"}))); + valid_params.set("Parameter Name","NOT SET BY USER"); + valid_params.set("Response Name","NOT SET BY USER"); + valid_params.set("Use FD Gradient",true); + + objective_params->validateParametersAndSetDefaults(valid_params); + + objective_type_ = objective_params->get("Objective Type"); + sensitivity_method_ = objective_params->get("Sensitivity Method"); + use_fd_gradient_ = objective_params->get("Use FD Gradient"); + + std::tie(model_,global_data_,mesh_,response_library_,stk_io_response_library_,lin_obj_factory_,global_indexer_) = + user_app::buildModelEvaluator(input_params_,comm_); + + param_index_ = std::get<0>(user_app::findParameterIndex(objective_params->get("Parameter Name"),*model_)); + response_index_ = std::get<0>(user_app::findResponseIndex(objective_params->get("Response Name"),*model_)); +} + +template +Real +TransientReducedObjective:: +value( const ROL::Vector &p, Real &tol ) +{ + using Teuchos::RCP; + typedef Thyra::ModelEvaluatorBase MEB; + + // Run tempus and compute response for specified parameter values + MEB::InArgs inArgs = model_->getNominalValues(); + MEB::OutArgs outArgs = model_->createOutArgs(); + const ROL::ThyraVector& thyra_p = + Teuchos::dyn_cast >(p); + inArgs.set_p(param_index_, thyra_p.getVector()); + RCP > g = + Thyra::createMember(model_->get_g_space(response_index_)); + outArgs.set_g(response_index_, g); + run_tempus(inArgs, outArgs); + + // Compute objective + if (target_ != Teuchos::null) + Thyra::Vp_StV(g.ptr(), -1.0, *target_); // g <- g - target + Real val = 0.0; + if (objective_type_ == "Sum of Squares") { + Thyra::ele_wise_scale(*g, g.ptr()); // g <- g.*g + val = Thyra::sum(*g); // val = \sum_i g_i + } + else if (objective_type_ == "Sum") + val = Thyra::sum(*g); // val = \sum_i g_i + else + TEUCHOS_TEST_FOR_EXCEPTION( + true, std::logic_error, + "Invalid objective type " << objective_type_ << "!" << std::endl + << "Valid choices are: \"Sum\" and \"Sum of Squares\"."); + + if (print_debug_) { + // const ROL::ThyraVector& thyra_p = Teuchos::dyn_cast >(*p); + // ROL::ThyraVector& thyra_g = Teuchos::dyn_cast >(*g); + *os_ << "TransientReducedObjective::value(): p = " << Thyra::get_ele(*(thyra_p.getVector()),0) + << ", g = " << Thyra::get_ele(*g,0) + << ", objective_val=" << val + << std::endl; + } + + return val; +} + +template +void +TransientReducedObjective:: +gradient(ROL::Vector &grad, const ROL::Vector &p, Real &tol) +{ + if (use_fd_gradient_) { + ROL::Objective::gradient(grad, p, tol); + return; + } + + using Teuchos::RCP; + typedef Thyra::ModelEvaluatorBase MEB; + + // Run tempus and compute response gradient for specified parameter values + const int num_p = model_->get_p_space(param_index_)->dim(); + const int num_g = model_->get_g_space(response_index_)->dim(); + MEB::InArgs inArgs = model_->getNominalValues(); + MEB::OutArgs outArgs = model_->createOutArgs(); + const ROL::ThyraVector& thyra_p = + Teuchos::dyn_cast >(p); + inArgs.set_p(param_index_, thyra_p.getVector()); + RCP > g = + Thyra::createMember(model_->get_g_space(response_index_)); + RCP > dgdp; + MEB::EDerivativeMultiVectorOrientation dgdp_orientation = + MEB::DERIV_MV_JACOBIAN_FORM; + if (sensitivity_method_ == "Adjoint"|| + sensitivity_method_ == "Pseudotransient Adjoint") { + // Adjoint, PseudoTransientAdjoint integrators require gradient form + // (but does not necessarily require the model to support gradient form) + dgdp = + Thyra::createMembers(model_->get_p_space(param_index_), num_g); + dgdp_orientation = MEB::DERIV_MV_GRADIENT_FORM; + } + else { + // Form determined by what model supports + MEB::DerivativeSupport support = + outArgs.supports(MEB::OUT_ARG_DgDp, response_index_, param_index_); + if (support.supports(MEB::DERIV_MV_JACOBIAN_FORM)) { + dgdp = + Thyra::createMembers(model_->get_g_space(response_index_), num_p); + dgdp_orientation = MEB::DERIV_MV_JACOBIAN_FORM; + } + else if (support.supports(MEB::DERIV_MV_GRADIENT_FORM)) { + dgdp = + Thyra::createMembers(model_->get_p_space(param_index_), num_g); + dgdp_orientation = MEB::DERIV_MV_GRADIENT_FORM; + } + else + TEUCHOS_TEST_FOR_EXCEPTION( + true, std::logic_error, + "Model must support Jacobian or gradient forms for dg/dp!"); + } + outArgs.set_DgDp(response_index_, param_index_, + MEB::DerivativeMultiVector(dgdp, dgdp_orientation)); + outArgs.set_g(response_index_, g); + run_tempus(inArgs, outArgs); + + // Compute objective gradient + RCP > final_grad = + Teuchos::dyn_cast >(grad).getVector()->col(0); + if (target_ != Teuchos::null) + Thyra::Vp_StV(g.ptr(), -1.0, *target_); // g <- g - target + if (dgdp_orientation == MEB::DERIV_MV_JACOBIAN_FORM) { + for (int j=0; j > dgdp_j = dgdp->col(j); + Real grad_val = 0.0; + if (objective_type_ == "Sum of Squares") { + Thyra::ele_wise_prod_update(2.0, *g, dgdp_j.ptr()); + grad_val = Thyra::sum(*dgdp_j); + } + else if (objective_type_ == "Sum") + grad_val = Thyra::sum(*dgdp_j); + else + TEUCHOS_TEST_FOR_EXCEPTION( + true, std::logic_error, + "Invalid objective type " << objective_type_ << "!" << std::endl + << "Valid choices are: \"Sum\" and \"Sum of Squares\"."); + Thyra::set_ele(j, grad_val, final_grad.ptr()); + } + } + else { // gradient form + Thyra::assign(final_grad.ptr(), 0.0); + for (int i=0; icol(i))); + } + else if (objective_type_ == "Sum") { + Thyra::Vp_StV(final_grad.ptr(), 1.0, *(dgdp->col(i))); + } + else + TEUCHOS_TEST_FOR_EXCEPTION( + true, std::logic_error, + "Invalid objective type " << objective_type_ << "!" << std::endl + << "Valid choices are: \"Sum\" and \"Sum of Squares\"."); + } + } +} + +template +void +TransientReducedObjective:: +set_target(const Teuchos::RCP >& target) { + target_ = target; +} + +template +void +TransientReducedObjective:: +set_target(const Teuchos::RCP >& target) { + using Teuchos::rcp_dynamic_cast; + target_ = + rcp_dynamic_cast >(target, true)->getVector(); +} + +template +Teuchos::RCP > +TransientReducedObjective:: +create_design_vector() const { + typedef Thyra::ModelEvaluatorBase MEB; + + Teuchos::RCP > p = + Thyra::createMember(model_->get_p_space(param_index_)); + MEB::InArgs nominalValues = model_->getNominalValues(); + if (nominalValues.get_p(param_index_) != Teuchos::null) + Thyra::assign(p.ptr(), *(nominalValues.get_p(param_index_))); + else + Thyra::assign(p.ptr(), Teuchos::ScalarTraits::zero()); + return Teuchos::rcp(new ROL::ThyraVector(p)); +} + +template +Teuchos::RCP > +TransientReducedObjective:: +create_response_vector() const { + Teuchos::RCP > g = + Thyra::createMember(model_->get_g_space(response_index_)); + Thyra::assign(g.ptr(), Teuchos::ScalarTraits::zero()); + return Teuchos::rcp(new ROL::ThyraVector(g)); +} + +template +void +TransientReducedObjective:: +run_tempus(ROL::Vector& r, const ROL::Vector& p) +{ + typedef Thyra::ModelEvaluatorBase MEB; + + MEB::InArgs inArgs = model_->getNominalValues(); + MEB::OutArgs outArgs = model_->createOutArgs(); + const ROL::ThyraVector& thyra_p = + Teuchos::dyn_cast >(p); + ROL::ThyraVector& thyra_r = + Teuchos::dyn_cast >(r); + inArgs.set_p(param_index_, thyra_p.getVector()); + outArgs.set_g(response_index_, thyra_r.getVector()); + run_tempus(inArgs, outArgs); +} + +template +void +TransientReducedObjective:: +run_tempus(const Thyra::ModelEvaluatorBase::InArgs& inArgs, + const Thyra::ModelEvaluatorBase::OutArgs& outArgs) +{ + using Teuchos::rcp; + using Teuchos::RCP; + using Teuchos::rcpFromRef; + typedef Thyra::ModelEvaluatorBase MEB; + typedef Thyra::DefaultNominalBoundsOverrideModelEvaluator DNBOME; + + // Override nominal values in model to supplied inArgs + RCP wrapped_model = rcp(new DNBOME(model_, rcpFromRef(inArgs))); + + Real t; + RCP > x, x_dot; + RCP > dxdp, dxdotdp; + RCP > g = outArgs.get_g(response_index_); + RCP > dgdp; + MEB::EDerivativeMultiVectorOrientation dgdp_orientation = Thyra::ModelEvaluatorBase::DERIV_MV_JACOBIAN_FORM; + bool dgdp_supported = outArgs.supports(Thyra::ModelEvaluatorBase::OUT_ARG_DgDp,response_index_,param_index_).supports(Thyra::ModelEvaluatorBase::DERIV_MV_JACOBIAN_FORM) || outArgs.supports(Thyra::ModelEvaluatorBase::OUT_ARG_DgDp,response_index_,param_index_).supports(Thyra::ModelEvaluatorBase::DERIV_MV_GRADIENT_FORM); + if (dgdp_supported) { + dgdp = outArgs.get_DgDp(response_index_, param_index_).getMultiVector(); + dgdp_orientation = outArgs.get_DgDp(response_index_, param_index_).getMultiVectorOrientation(); + } + + // Create and run integrator + if (dgdp != Teuchos::null && sensitivity_method_ == "Forward") { + RCP > integrator = + Tempus::createIntegratorForwardSensitivity(tempus_params_, wrapped_model); + const bool integratorStatus = integrator->advanceTime(); + TEUCHOS_TEST_FOR_EXCEPTION( + !integratorStatus, std::logic_error, "Integrator failed!"); + + // Get final state + t = integrator->getTime(); + x = integrator->getX(); + x_dot = integrator->getXDot(); + dxdp = integrator->getDxDp(); + dxdotdp = integrator->getDXDotDp(); + } + else if (dgdp != Teuchos::null && sensitivity_method_ == "Adjoint") { + RCP > integrator = + Tempus::createIntegratorAdjointSensitivity(tempus_params_, wrapped_model); + const bool integratorStatus = integrator->advanceTime(); + TEUCHOS_TEST_FOR_EXCEPTION( + !integratorStatus, std::logic_error, "Integrator failed!"); + + // Get final state + t = integrator->getTime(); + x = integrator->getX(); + x_dot = integrator->getXDot(); + Thyra::assign(dgdp.ptr(), *(integrator->getDgDp())); + } + else if (dgdp != Teuchos::null && + sensitivity_method_ == "Pseudotransient Forward") { + RCP > integrator = + Tempus::createIntegratorPseudoTransientForwardSensitivity(tempus_params_, + wrapped_model); + const bool integratorStatus = integrator->advanceTime(); + TEUCHOS_TEST_FOR_EXCEPTION( + !integratorStatus, std::logic_error, "Integrator failed!"); + + // Get final state + t = integrator->getTime(); + x = integrator->getX(); + x_dot = integrator->getXDot(); + dxdp = integrator->getDxDp(); + dxdotdp = integrator->getDXDotDp(); + } + else if (dgdp != Teuchos::null && + sensitivity_method_ == "Pseudotransient Adjoint") { + RCP > integrator = + Tempus::integratorPseudoTransientAdjointSensitivity(tempus_params_, + wrapped_model); + const bool integratorStatus = integrator->advanceTime(); + TEUCHOS_TEST_FOR_EXCEPTION( + !integratorStatus, std::logic_error, "Integrator failed!"); + + // Get final state + t = integrator->getTime(); + x = integrator->getX(); + x_dot = integrator->getXDot(); + Thyra::assign(dgdp.ptr(), *(integrator->getDgDp())); + } + else if (dgdp != Teuchos::null) { + TEUCHOS_TEST_FOR_EXCEPTION( + true, std::logic_error, + "Invalid sensitivity method " << sensitivity_method_ << "!" << std::endl + << "Valid choices are: Forward, Adjoint, Pseudotransient Forward, " + << " and Pseudotransient Adjoint."); + } + else { + /* + RCP > integrator = + Tempus::createIntegratorBasic(tempus_params_, wrapped_model); + const bool integratorStatus = integrator->advanceTime(); + TEUCHOS_TEST_FOR_EXCEPTION( + !integratorStatus, std::logic_error, "Integrator failed!"); + */ + + auto input_params_copy = Teuchos::parameterList(*input_params_); + const bool override_nox_output = true; + auto integrator = user_app::buildTimeIntegrator(input_params_copy,comm_,wrapped_model,mesh_,response_library_, + stk_io_response_library_, + lin_obj_factory_,global_indexer_, + override_nox_output); + + bool integratorStatus = integrator->advanceTime(); + TEUCHOS_ASSERT(integratorStatus); + + // Get final state + t = integrator->getTime(); + x = integrator->getX(); + x_dot = integrator->getXDot(); + } + + // Evaluate response at final state + const int num_g = model_->get_g_space(response_index_)->dim(); + MEB::InArgs modelInArgs = inArgs; + MEB::OutArgs modelOutArgs = outArgs; + modelInArgs.set_x(x); + if (modelInArgs.supports(MEB::IN_ARG_x_dot)) modelInArgs.set_x_dot(x_dot); + if (modelInArgs.supports(MEB::IN_ARG_t)) modelInArgs.set_t(t); + RCP > dgdx, dgdxdot; + MEB::EDerivativeMultiVectorOrientation dgdx_orientation = + MEB::DERIV_MV_JACOBIAN_FORM; + MEB::EDerivativeMultiVectorOrientation dgdxdot_orientation = + MEB::DERIV_MV_JACOBIAN_FORM; + if (dgdp != Teuchos::null && + (sensitivity_method_ == "Forward" || + sensitivity_method_ == "Pseudotransient Forward")) { + MEB::DerivativeSupport dgdx_support = + outArgs.supports(MEB::OUT_ARG_DgDx, response_index_); + if (!dgdx_support.none()) { + if (dgdx_support.supports(MEB::DERIV_MV_GRADIENT_FORM)) { + dgdx = Thyra::createMembers(model_->get_x_space(), num_g); + dgdx_orientation = MEB::DERIV_MV_GRADIENT_FORM; + } + else + TEUCHOS_TEST_FOR_EXCEPTION( + true, std::logic_error, + "Model must support gradient forms for dg/dx!"); + modelOutArgs.set_DgDx( + response_index_, + MEB::DerivativeMultiVector(dgdx, dgdx_orientation)); + } + + MEB::DerivativeSupport dgdxdot_support = + modelOutArgs.supports(MEB::OUT_ARG_DgDx_dot, response_index_); + if (!dgdxdot_support.none()) { + if (dgdxdot_support.supports(MEB::DERIV_MV_GRADIENT_FORM)) { + dgdxdot = Thyra::createMembers(model_->get_x_space(), num_g); + dgdxdot_orientation = MEB::DERIV_MV_GRADIENT_FORM; + } + else + TEUCHOS_TEST_FOR_EXCEPTION( + true, std::logic_error, + "Model must support Jacobian or gradient forms for dg/dx!"); + modelOutArgs.set_DgDx_dot( + response_index_, + MEB::DerivativeMultiVector(dgdxdot, dgdxdot_orientation)); + } + } + else if (dgdp != Teuchos::null && + (sensitivity_method_ == "Adjoint" || + sensitivity_method_ == "Pseudotransient Adjoint")) { + // Clear dg/dp as an out arg since it was already computed by the adjoint + // integrator + modelOutArgs.set_DgDp(response_index_, param_index_, + MEB::Derivative()); + } + + model_->evalModel(modelInArgs, modelOutArgs); + + // dg/dp = dg/dp + dg/dx*dx/dp + dg/dx_dot*dx_dot/dp + // We assume dg/dx, dg/dxdot are in gradient form while dxdp, dxdotdp are in + // Jacobian form + if (dgdp != Teuchos::null && dgdx != Teuchos::null) { + if (dgdp_orientation == MEB::DERIV_MV_JACOBIAN_FORM) + dgdx->apply(Thyra::TRANS, *dxdp, dgdp.ptr(), Real(1.0), Real(1.0)); + else + dxdp->apply(Thyra::TRANS, *dgdx, dgdp.ptr(), Real(1.0), Real(1.0)); + } + if (dgdp != Teuchos::null && dgdxdot != Teuchos::null) { + if (dgdp_orientation == MEB::DERIV_MV_JACOBIAN_FORM) + dgdxdot->apply(Thyra::TRANS, *dxdotdp, dgdp.ptr(), Real(1.0), Real(1.0)); + else + dxdotdp->apply(Thyra::TRANS, *dgdxdot, dgdp.ptr(), Real(1.0), Real(1.0)); + } +} + +} // namespace ROL diff --git a/packages/panzer/adapters-stk/example/main_driver/user_app_Utilities.cpp b/packages/panzer/adapters-stk/example/main_driver/user_app_Utilities.cpp new file mode 100644 index 000000000000..2d4dadc65a78 --- /dev/null +++ b/packages/panzer/adapters-stk/example/main_driver/user_app_Utilities.cpp @@ -0,0 +1,309 @@ +// @HEADER +// ***************************************************************************** +// Panzer: A partial differential equation assembly +// engine for strongly coupled complex multiphysics systems +// +// Copyright 2011 NTESS and the Panzer contributors. +// SPDX-License-Identifier: BSD-3-Clause +// ***************************************************************************** +// @HEADER + +#include "Teuchos_ConfigDefs.hpp" +#include "Teuchos_RCP.hpp" +#include "Teuchos_TimeMonitor.hpp" +#include "Teuchos_StackedTimer.hpp" +#include "Teuchos_DefaultComm.hpp" +#include "Teuchos_CommHelpers.hpp" +#include "Teuchos_GlobalMPISession.hpp" +#include "Teuchos_CommandLineProcessor.hpp" +#include "Teuchos_XMLParameterListHelpers.hpp" +#include "Teuchos_YamlParameterListHelpers.hpp" +#include "Teuchos_FancyOStream.hpp" +#include "Teuchos_oblackholestream.hpp" +#include "Teuchos_Assert.hpp" +#include "Teuchos_as.hpp" + +#include "Panzer_NodeType.hpp" +#include "PanzerAdaptersSTK_config.hpp" +#include "Panzer_STK_ModelEvaluatorFactory.hpp" +#include "Panzer_ClosureModel_Factory_TemplateManager.hpp" +#include "Panzer_String_Utilities.hpp" +#include "Panzer_ThyraObjContainer.hpp" +#include "Thyra_VectorSpaceBase.hpp" + +#include "NOX_Utils.H" +#include "NOX_Observer_Print.hpp" + +#include "Tempus_IntegratorBasic.hpp" +#include "Tempus_IntegratorObserverComposite.hpp" +#include "Tempus_IntegratorObserverBasic.hpp" +#include "Tempus_StepperFactory.hpp" + +#include "user_app_ClosureModel_Factory_TemplateBuilder.hpp" +#include "user_app_EquationSetFactory.hpp" +#include "user_app_BCStrategy_Factory.hpp" +#include "user_app_NOXObserverFactory.hpp" +#ifdef PANZER_HAVE_TEMPUS +#include "user_app_TempusObserverFactory.hpp" +#endif +#include "user_app_ResponseEvaluatorFactory_HOFlux.hpp" + +#include "Panzer_ROLTempusReducedObjective.hpp" +#include "user_app_TransientObjective.hpp" + +#include "Panzer_ResponseEvaluatorFactory_Probe.hpp" + +#include "ROL_OptimizationProblem.hpp" +#include "ROL_OptimizationSolver.hpp" +#include "ROL_RandomVector.hpp" + +#include + +#include +#include +#include + +#include "user_app_Utilities.hpp" + +void user_app::addResponsesToModelEvaluatorFactory(const Teuchos::ParameterList& responses, + panzer_stk::ModelEvaluatorFactory& me_factory) +{ + for(Teuchos::ParameterList::ConstIterator itr=responses.begin();itr!=responses.end();++itr) { + const std::string name = responses.name(itr); + TEUCHOS_ASSERT(responses.entry(itr).isList()); + Teuchos::ParameterList & lst = Teuchos::getValue(responses.entry(itr)); + + if (lst.get("Type") == "Functional") { + + panzer::FunctionalResponse_Builder builder; + builder.comm = MPI_COMM_WORLD; + builder.cubatureDegree = 2; + builder.requiresCellIntegral = true; + builder.quadPointField = lst.get("Field Name"); + + std::vector eblocks; + panzer::StringTokenizer(eblocks,lst.get("Element Blocks"),",",true); + + std::vector wkst_descs; + for(std::size_t i=0;i("Type") == "Point Value") { + panzer::ProbeResponse_Builder builder; + builder.comm = MPI_COMM_WORLD; + builder.point = Teuchos::Array{0.5,0.5}; + builder.cubatureDegree = 2; + builder.fieldName = lst.get("Field Name"); + builder.applyDirichletToDerivative = false; + + std::vector eblocks; + panzer::StringTokenizer(eblocks,lst.get("Element Blocks"),",",true); + + std::vector descriptors; + for(std::size_t i=0;i("Type") << "\" is not supported!"); + } + } +} + +std::tuple< Teuchos::RCP>, + Teuchos::RCP, + Teuchos::RCP, + Teuchos::RCP>, + Teuchos::RCP>, + Teuchos::RCP>, + Teuchos::RCP + > +user_app::buildModelEvaluator(const Teuchos::RCP& input_params, + const Teuchos::RCP>& comm) +{ + using namespace Teuchos; + using namespace Tempus; + + Teuchos::ParameterList solver_factories = input_params->sublist("Solver Factories"); + input_params->remove("Solver Factories"); + + // Add in the application specific equation set factory + Teuchos::RCP eqset_factory = Teuchos::rcp(new user_app::MyFactory); + + // Add in the application specific closure model factory + user_app::MyModelFactory_TemplateBuilder cm_builder; + panzer::ClosureModelFactory_TemplateManager cm_factory; + cm_factory.buildObjects(cm_builder); + + // Add in the application specific bc factory + user_app::BCFactory bc_factory; + + // Create the global data + Teuchos::RCP global_data = panzer::createGlobalData(); + input_params->sublist("User Data").set("Comm", comm); // Required for GlobalStatistics + + Teuchos::RCP > stkIOResponseLibrary + = Teuchos::rcp(new panzer::ResponseLibrary()); + + Teuchos::RCP> physics; + Teuchos::RCP> solver; + Teuchos::RCP> rLibrary; + std::vector> physicsBlocks; + Teuchos::RCP> linObjFactory; + Teuchos::RCP globalIndexer; + Teuchos::RCP mesh; + { + Teuchos::ParameterList responses = input_params->sublist("Responses"); + input_params->remove("Responses"); + + panzer_stk::ModelEvaluatorFactory me_factory; + me_factory.setParameterList(input_params); + auto strat_list = Teuchos::parameterList(); + *strat_list = input_params->sublist("Solution Control",true).sublist("Tempus",true).sublist("My Example Stepper",true).sublist("My Example Solver",true).sublist("NOX",true).sublist("Direction",true).sublist("Newton",true).sublist("Stratimikos Linear Solver",true).sublist("Stratimikos",true); + me_factory.setStratimikosList(strat_list); + me_factory.buildObjects(comm,global_data,eqset_factory,bc_factory,cm_factory); + + user_app::addResponsesToModelEvaluatorFactory(responses,me_factory); + me_factory.buildResponses(cm_factory); + + physicsBlocks = me_factory.getPhysicsBlocks(); + physics = me_factory.getPhysicsModelEvaluator(); + rLibrary = me_factory.getResponseLibrary(); + linObjFactory = me_factory.getLinearObjFactory(); + globalIndexer = me_factory.getGlobalIndexer(); + mesh = me_factory.getMesh(); + } + + // setup outputs to mesh on the stkIOResponseLibrary + { + stkIOResponseLibrary->initialize(*rLibrary); + + Teuchos::ParameterList user_data(input_params->sublist("User Data")); + user_data.set("Workset Size",input_params->sublist("Assembly").get("Workset Size")); + + std::vector eBlocks; + mesh->getElementBlockNames(eBlocks); + + panzer_stk::RespFactorySolnWriter_Builder builder; + builder.mesh = mesh; + stkIOResponseLibrary->addResponse("Main Field Output",eBlocks,builder); + + stkIOResponseLibrary->buildResponseEvaluators(physicsBlocks, + cm_factory, + input_params->sublist("Closure Models"), + user_data); + } + + return {physics,global_data,mesh,rLibrary,stkIOResponseLibrary,linObjFactory,globalIndexer}; +} + +Teuchos::RCP> +user_app::buildTimeIntegrator(const Teuchos::RCP& input_params, + const Teuchos::RCP>& comm, + Teuchos::RCP> physics, + Teuchos::RCP mesh, + Teuchos::RCP> rLibrary, + Teuchos::RCP> stkIOResponseLibrary, + Teuchos::RCP> linObjFactory, + Teuchos::RCP globalIndexer, + const bool overrideNoxOutput) +{ + using namespace Teuchos; + using namespace Tempus; + + RCP tempus_pl = parameterList(input_params->sublist("Solution Control").sublist("Tempus")); + + if (overrideNoxOutput) { + auto nox_utils = Teuchos::rcp(new NOX::Utils(NOX::Utils::Error,comm->getRank(),0,3)); + auto print_nonlinear = Teuchos::rcp(new NOX::ObserverPrint(nox_utils,10)); + tempus_pl->sublist("My Example Stepper",true) + .sublist("My Example Solver",true) + .sublist("NOX",true) + .sublist("Solver Options",true) + .set>("Observer",print_nonlinear); + } + + // Tempus is overriding NOX parameters with its own defaults. Need + // to fix this. It also does not validate because of arbitrary + // solver name. Should validate but use + // disableRecursiveValidation() to avoid errors with arbitrary + // names. + + const bool doInitialization = false; + auto integrator = createIntegratorBasic(tempus_pl, physics, doInitialization); + + RCP noxList = parameterList("Correct NOX Params"); + *noxList = tempus_pl->sublist("My Example Stepper",true).sublist("My Example Solver",true).sublist("NOX",true); + // noxList->print(std::cout); + integrator->getStepper()->getSolver()->setParameterList(noxList); + integrator->initialize(); + + // Setting observers on tempus breaks the screen output of the + // time steps! It replaces IntegratorObserverBasic which handles + // IO. Is this at least documented? + { + RCP> tempus_observers = rcp(new Tempus::IntegratorObserverComposite); + RCP tof = + Teuchos::rcp(new user_app::TempusObserverFactory(stkIOResponseLibrary,rLibrary->getWorksetContainer())); + tempus_observers->addObserver(Teuchos::rcp(new Tempus::IntegratorObserverBasic)); + tempus_observers->addObserver(tof->buildTempusObserver(mesh,globalIndexer,linObjFactory)); + integrator->setObserver(tempus_observers); + } + + RCP> x0 = physics->getNominalValues().get_x()->clone_v(); + integrator->initializeSolutionHistory(0.0, x0); + + return integrator; +} + +std::tuple user_app::findParameterIndex(const std::string& p_name,const Thyra::ModelEvaluator& me) +{ + bool found = false; + int index = -1; + int sub_index = -1; + for (int p = 0; p < me.Np(); ++p) { + auto p_names = me.get_p_names(p); + for (Teuchos::Ordinal p_sub=0; p_sub < p_names->size(); ++p_sub) { + if ((*p_names)[p_sub] == p_name) { + found = true; + index = p; + sub_index = p_sub; + break; + } + } + if (found) break; + } + + TEUCHOS_TEST_FOR_EXCEPTION(!found,std::runtime_error, + "ERROR: the parameter \"" << p_name << "\" is not a valid parameter name in the model evaluator!"); + + return {index,sub_index}; +} + +std::tuple user_app::findResponseIndex(const std::string& g_name,const Thyra::ModelEvaluator& me) +{ + bool found = false; + int index = -1; + int sub_index = -1; + for (int g = 0; g < me.Ng(); ++g) { + auto g_names = me.get_g_names(g); + for (Teuchos::Ordinal g_sub=0; g_sub < g_names.size(); ++g_sub) { + std::cout << "g(" << g << "," << g_sub << ")=" << g_names[g_sub] << std::endl; + if (g_names[g_sub] == g_name) { + found = true; + index = g; + sub_index = g_sub; + break; + } + } + if (found) break; + } + + TEUCHOS_TEST_FOR_EXCEPTION(!found,std::runtime_error, + "ERROR: the response \"" << g_name << "\" is not a valid response name in the model evaluator!"); + + return {index,sub_index}; +} diff --git a/packages/panzer/adapters-stk/example/main_driver/user_app_Utilities.hpp b/packages/panzer/adapters-stk/example/main_driver/user_app_Utilities.hpp new file mode 100644 index 000000000000..14aeca4a2494 --- /dev/null +++ b/packages/panzer/adapters-stk/example/main_driver/user_app_Utilities.hpp @@ -0,0 +1,145 @@ +// @HEADER +// ***************************************************************************** +// Panzer: A partial differential equation assembly +// engine for strongly coupled complex multiphysics systems +// +// Copyright 2011 NTESS and the Panzer contributors. +// SPDX-License-Identifier: BSD-3-Clause +// ***************************************************************************** +// @HEADER + +#ifndef PANZER_ADAPTERS_STK_MAIN_DRIVER_USER_APP_UTILITIES_HPP +#define PANZER_ADAPTERS_STK_MAIN_DRIVER_USER_APP_UTILITIES_HPP + +#include "Teuchos_ConfigDefs.hpp" +#include "Teuchos_RCP.hpp" +#include "Teuchos_DefaultComm.hpp" +#include "Teuchos_CommHelpers.hpp" + +#include "Panzer_NodeType.hpp" +#include "PanzerAdaptersSTK_config.hpp" +#include "Panzer_STK_ModelEvaluatorFactory.hpp" +#include "Panzer_ClosureModel_Factory_TemplateManager.hpp" +#include "Panzer_String_Utilities.hpp" +#include "Panzer_ThyraObjContainer.hpp" +#include "Thyra_VectorSpaceBase.hpp" + +#include "NOX_Utils.H" +#include "NOX_Observer_Print.hpp" + +#include "Tempus_IntegratorBasic.hpp" + +#include "user_app_ClosureModel_Factory_TemplateBuilder.hpp" +#include "user_app_EquationSetFactory.hpp" +#include "user_app_BCStrategy_Factory.hpp" +#include "user_app_NOXObserverFactory.hpp" +#include "user_app_TempusObserverFactory.hpp" +#include "user_app_ResponseEvaluatorFactory_HOFlux.hpp" + +#include +#include + +namespace user_app { + + void addResponsesToModelEvaluatorFactory(const Teuchos::ParameterList& response_sublist, + panzer_stk::ModelEvaluatorFactory& me_factory); + + void computeAndPrintResponses(const Teuchos::RCP>& me, + const auto& x, + const auto& x_dot, + const auto& t, + const auto& global_data, + std::ostream& out); + + std::tuple< Teuchos::RCP>, + Teuchos::RCP, + Teuchos::RCP, + Teuchos::RCP>, // normal response lib + Teuchos::RCP>, // stk io response lib + Teuchos::RCP>, + Teuchos::RCP + > + buildModelEvaluator(const Teuchos::RCP& input_params, + const Teuchos::RCP>& comm); + + Teuchos::RCP> + buildTimeIntegrator(const Teuchos::RCP& input_params, + const Teuchos::RCP>& comm, + Teuchos::RCP> me, + Teuchos::RCP mesh, + Teuchos::RCP> rLibrary, + Teuchos::RCP> stkIOResponseLibrary, + Teuchos::RCP> linObjFactory, + Teuchos::RCP globalIndexer, + const bool overrideNoxOutput); + + std::tuple findParameterIndex(const std::string& p_name,const Thyra::ModelEvaluator& me); + + std::tuple findResponseIndex(const std::string& g_name,const Thyra::ModelEvaluator& me); +} + +void user_app::computeAndPrintResponses(const Teuchos::RCP>& physics, + const auto& x, + const auto& x_dot, + const auto& t, + const auto& global_data, + std::ostream& os) +{ + if(physics->Ng()>0) { + + os << "Ng = " << physics->Ng() << std::endl; + os << "Np = " << physics->Np() << std::endl; + Thyra::ModelEvaluatorBase::InArgs respInArgs = physics->createInArgs(); + Thyra::ModelEvaluatorBase::OutArgs respOutArgs = physics->createOutArgs(); + + TEUCHOS_ASSERT(physics->Ng()==respOutArgs.Ng()); + + respInArgs.set_x(x); + respInArgs.set_x_dot(x_dot); + + for(int g=0;g > response = Thyra::createMember(*physics->get_g_space(g)); + respOutArgs.set_g(g,response); + + for (int p = 0; p < physics->Np(); ++p) { + bool derivative_supported = respOutArgs.supports(Thyra::ModelEvaluatorBase::OUT_ARG_DgDp,g,p).supports(Thyra::ModelEvaluatorBase::DERIV_MV_JACOBIAN_FORM); + os << "DgDp(" << g << "," << p << ") supports = " << derivative_supported << std::endl; + if (derivative_supported) { + auto dgdp = physics->create_DgDp_op(g,p); + respOutArgs.set_DgDp(g,p,Thyra::ModelEvaluatorBase::Derivative(dgdp)); + } + } + } + + physics->evalModel(respInArgs, respOutArgs); + + { + auto parameter_library = global_data->pl; + TEUCHOS_ASSERT(nonnull(parameter_library)); + for (auto i=parameter_library->begin();i != parameter_library->end(); ++i) { + os << "Sacado::ParameterLibrary: " << i->first << std::endl; + } + } + + { + auto nominalValues = physics->getNominalValues(); + for (int i=0; i < respInArgs.Np();++i) { + auto p = nominalValues.get_p(i); + auto p_names = physics->get_p_names(i); + os << "ModelEvaluator::NominalValues Parameter Value: \"" << (*p_names)[0] << "\" = " << Thyra::get_ele(*p,0) << std::endl; + } + } + + for (int i=0;i > response = respOutArgs.get_g(i); + TEUCHOS_ASSERT(response!=Teuchos::null); + os << "Response Value: \"" << physics->get_g_names(i)[0] << "\" = " << Thyra::get_ele(*response,0) << std::endl; + for (int j=0; j < respOutArgs.Np(); ++j) { + // os << " dg(" << i << ")/dp(" << j << ") supports(GRAD_FORM) = " << respOutArgs.supports(Thyra::ModelEvaluatorBase::OUT_ARG_DgDp,i,j).supports(Thyra::ModelEvaluatorBase::DERIV_MV_GRADIENT_FORM) << std::endl; + // os << " dg(" << i << ")/dp(" << j << ") supports(JAC_FORM) = " << respOutArgs.supports(Thyra::ModelEvaluatorBase::OUT_ARG_DgDp,i,j).supports(Thyra::ModelEvaluatorBase::DERIV_MV_JACOBIAN_FORM) << std::endl; + } + } + } +} + +#endif diff --git a/packages/panzer/disc-fe/src/Panzer_ModelEvaluator_impl.hpp b/packages/panzer/disc-fe/src/Panzer_ModelEvaluator_impl.hpp index 1f3101881359..92d60c06fcc1 100644 --- a/packages/panzer/disc-fe/src/Panzer_ModelEvaluator_impl.hpp +++ b/packages/panzer/disc-fe/src/Panzer_ModelEvaluator_impl.hpp @@ -1647,6 +1647,7 @@ void panzer::ModelEvaluator:: evalModelImpl_basic_g(const Thyra::ModelEvaluatorBase::InArgs &inArgs, const Thyra::ModelEvaluatorBase::OutArgs &outArgs) const { + PANZER_FUNC_TIME_MONITOR("panzer::ModelEvaluator::evalModelImpl_basic_g()"); // optional sanity check // TEUCHOS_ASSERT(required_basic_g(outArgs)); @@ -1684,6 +1685,7 @@ panzer::ModelEvaluator:: evalModelImpl_basic_dgdx(const Thyra::ModelEvaluatorBase::InArgs &inArgs, const Thyra::ModelEvaluatorBase::OutArgs &outArgs) const { + PANZER_FUNC_TIME_MONITOR("panzer::ModelEvaluator::evalModelImpl_basic_dgdx()"); typedef Thyra::ModelEvaluatorBase MEB; // optional sanity check @@ -1730,6 +1732,7 @@ panzer::ModelEvaluator:: evalModelImpl_basic_dgdp_scalar(const Thyra::ModelEvaluatorBase::InArgs &inArgs, const Thyra::ModelEvaluatorBase::OutArgs &outArgs) const { + PANZER_FUNC_TIME_MONITOR("panzer::ModelEvaluator::evalModelImpl_basic_dgdp_scalar()"); using Teuchos::RCP; using Teuchos::rcp; using Teuchos::rcp_dynamic_cast; @@ -1818,6 +1821,7 @@ panzer::ModelEvaluator:: evalModelImpl_basic_dgdp_distro(const Thyra::ModelEvaluatorBase::InArgs &inArgs, const Thyra::ModelEvaluatorBase::OutArgs &outArgs) const { + PANZER_FUNC_TIME_MONITOR("panzer::ModelEvaluator::evalModelImpl_basic_dgdp_distro()"); typedef Thyra::ModelEvaluatorBase MEB; // optional sanity check @@ -1877,6 +1881,7 @@ panzer::ModelEvaluator:: evalModelImpl_basic_dfdp_scalar(const Thyra::ModelEvaluatorBase::InArgs &inArgs, const Thyra::ModelEvaluatorBase::OutArgs &outArgs) const { + PANZER_FUNC_TIME_MONITOR("panzer::ModelEvaluator::evalModelImpl_basic_dfdp_scalar()"); using Teuchos::RCP; using Teuchos::rcp_dynamic_cast; @@ -1982,7 +1987,7 @@ evalModelImpl_basic_dfdp_scalar(const Thyra::ModelEvaluatorBase::InArgs /////////////////////////////////////////////////////////////////////////////////////// if(totalParameterCount>0) { - PANZER_FUNC_TIME_MONITOR("panzer::ModelEvaluator::evalModel(df/dp)"); + PANZER_FUNC_TIME_MONITOR_DIFF("panzer::ModelEvaluator::evalModel(df/dp)",dfdp_eval); ae_tm_.getAsObject()->evaluate(ae_inargs); } } @@ -1993,7 +1998,7 @@ panzer::ModelEvaluator:: evalModelImpl_basic_dfdp_scalar_fd(const Thyra::ModelEvaluatorBase::InArgs &inArgs, const Thyra::ModelEvaluatorBase::OutArgs &outArgs) const { - PANZER_FUNC_TIME_MONITOR("panzer::ModelEvaluator::evalModel(df/dp)"); + PANZER_FUNC_TIME_MONITOR("panzer::ModelEvaluator::evalModelImpl_basic_dfdp_scalar_fd()"); using Teuchos::RCP; using Teuchos::rcp_dynamic_cast; @@ -2097,6 +2102,7 @@ panzer::ModelEvaluator:: evalModelImpl_basic_dfdp_distro(const Thyra::ModelEvaluatorBase::InArgs &inArgs, const Thyra::ModelEvaluatorBase::OutArgs &outArgs) const { + PANZER_FUNC_TIME_MONITOR("panzer::ModelEvaluator::evalModelImpl_basic_dfdp_distro()"); using Teuchos::RCP; using Teuchos::rcp_dynamic_cast; using Teuchos::null; diff --git a/packages/panzer/disc-fe/src/evaluators/Panzer_Parameter_impl.hpp b/packages/panzer/disc-fe/src/evaluators/Panzer_Parameter_impl.hpp index ce108d8ab80a..6c77f855f345 100644 --- a/packages/panzer/disc-fe/src/evaluators/Panzer_Parameter_impl.hpp +++ b/packages/panzer/disc-fe/src/evaluators/Panzer_Parameter_impl.hpp @@ -49,13 +49,12 @@ evaluateFields(typename TRAITS::EvalData workset) auto param_val = param->getValue(); auto target_field_v = target_field.get_static_view(); auto target_field_h = Kokkos::create_mirror_view(target_field_v); - + for (int cell=0; cell < workset.num_cells; ++cell) { for (std::size_t pt=0; pt