diff --git a/src/Albany_ObserverImpl.cpp b/src/Albany_ObserverImpl.cpp index 00356b6864..ea6a2253df 100644 --- a/src/Albany_ObserverImpl.cpp +++ b/src/Albany_ObserverImpl.cpp @@ -22,16 +22,12 @@ ObserverImpl (const Teuchos::RCP &app) void ObserverImpl:: observeSolution(double stamp, - const Thyra_Vector& nonOverlappedSolution, - const Teuchos::Ptr& nonOverlappedSolution_dxdp, - const Teuchos::Ptr& nonOverlappedSolutionDot, - const Teuchos::Ptr& nonOverlappedSolutionDotDot) + const Thyra_Vector& x, + const Teuchos::Ptr& dxdp, + const Teuchos::Ptr& x_dot, + const Teuchos::Ptr& x_dotdot) { - app_->evaluateStateFieldManager (stamp, - nonOverlappedSolution, - nonOverlappedSolutionDot, - nonOverlappedSolutionDotDot, - nonOverlappedSolution_dxdp); + app_->evaluateStateFieldManager (stamp, x, x_dot, x_dotdot, dxdp); //! update distributed parameters in the mesh auto distParamLib = app_->getDistributedParameterLibrary(); @@ -44,21 +40,28 @@ observeSolution(double stamp, } StatelessObserverImpl::observeSolution (stamp, - nonOverlappedSolution, - nonOverlappedSolution_dxdp, - nonOverlappedSolutionDot, - nonOverlappedSolutionDotDot); + Teuchos::rcpFromRef(x), + Teuchos::rcpFromPtr(x_dot), + Teuchos::rcpFromPtr(x_dotdot), + Teuchos::rcpFromPtr(dxdp)); } void ObserverImpl:: observeSolution(double stamp, - const Thyra_MultiVector& nonOverlappedSolution, - const Teuchos::Ptr& nonOverlappedSolution_dxdp) + const Thyra_MultiVector& x, + const Teuchos::Ptr& dxdp) { - app_->evaluateStateFieldManager(stamp, nonOverlappedSolution, - nonOverlappedSolution_dxdp); - StatelessObserverImpl::observeSolution(stamp, nonOverlappedSolution, - nonOverlappedSolution_dxdp); + app_->evaluateStateFieldManager(stamp, x, dxdp); + Teuchos::Ptr x_dot,x_dotdot; + + int x_ncols = x.domain()->dim(); + if (x_ncols==1) { + x_dot = x.col(1).ptr(); + } else if (x_ncols==2) { + x_dot = x.col(1).ptr(); + x_dotdot = x.col(2).ptr(); + } + observeSolution (stamp, *x.col(0), dxdp, x_dot, x_dotdot); } void ObserverImpl:: diff --git a/src/Albany_ObserverImpl.hpp b/src/Albany_ObserverImpl.hpp index f0a70c7609..399199c420 100644 --- a/src/Albany_ObserverImpl.hpp +++ b/src/Albany_ObserverImpl.hpp @@ -13,25 +13,26 @@ namespace Albany { -class ObserverImpl : public StatelessObserverImpl, public Piro::ROL_ObserverBase { +class ObserverImpl : public StatelessObserverImpl, + public Piro::ROL_ObserverBase +{ public: explicit ObserverImpl(const Teuchos::RCP& app); // Bring in scope the version we don't override using StatelessObserverImpl::observeSolution; - void observeSolution( - double stamp, const Thyra_Vector& nonOverlappedSolution, - const Teuchos::Ptr& nonOverlappedSolution_dxdp, - const Teuchos::Ptr& nonOverlappedSolutionDot, - const Teuchos::Ptr& nonOverlappedSolutionDotDot) override; + // Override from ROL_ObserverBase + void observeSolution (double stamp, const Thyra_Vector& x, + const Teuchos::Ptr& x_dxdp, + const Teuchos::Ptr& x_dot, + const Teuchos::Ptr& x_dotdot) override; - void observeSolution( - double stamp, const Thyra_MultiVector& nonOverlappedSolution, - const Teuchos::Ptr& nonOverlappedSolution_dxdp) override; + void observeSolution (double stamp, + const Thyra_MultiVector& x, + const Teuchos::Ptr& dxdp) override; - void parameterChanged( - const std::string& param) override; + void parameterChanged (const std::string& param) override; void parametersChanged() override; diff --git a/src/Albany_PiroObserver.cpp b/src/Albany_PiroObserver.cpp index a226323c4d..17370dbfb7 100644 --- a/src/Albany_PiroObserver.cpp +++ b/src/Albany_PiroObserver.cpp @@ -38,334 +38,222 @@ PiroObserver(const Teuchos::RCP &app, } void PiroObserver:: -observeSolution(const Thyra_Vector &solution) -{ - this->observeSolutionImpl(solution, Teuchos::ScalarTraits::zero()); - stepper_counter_++; -} - -void PiroObserver:: -observeSolution(const Thyra_Vector &solution, - const Thyra_MultiVector &solution_dxdp) -{ - this->observeSolutionImpl(solution, solution_dxdp, Teuchos::ScalarTraits::zero()); - stepper_counter_++; -} - -void PiroObserver:: -observeSolution(const Thyra_Vector &solution, +observeSolution(const Thyra_Vector& x, const ST stamp) { - this->observeSolutionImpl(solution, stamp); - stepper_counter_++; + this->observeSolutionImpl(Teuchos::rcpFromRef(x), + Teuchos::null, // xdot + Teuchos::null, // xdotdot + Teuchos::null, // dxdp + stamp); } void PiroObserver:: -observeSolution(const Thyra_Vector &solution, - const Thyra_MultiVector &solution_dxdp, +observeSolution(const Thyra_Vector &x, + const Thyra_MultiVector& dxdp, const ST stamp) { - this->observeSolutionImpl(solution, solution_dxdp, stamp); - stepper_counter_++; + this->observeSolutionImpl(Teuchos::rcpFromRef(x), + Teuchos::null, // xdot + Teuchos::null, // xdotdot + Teuchos::rcpFromRef(dxdp), + stamp); } void PiroObserver:: -observeSolution(const Thyra_Vector &solution, - const Thyra_Vector &solution_dot, +observeSolution(const Thyra_Vector& x, + const Thyra_Vector& x_dot, const ST stamp) { - this->observeSolutionImpl(solution, solution_dot, stamp); - stepper_counter_++; + this->observeSolutionImpl(Teuchos::rcpFromRef(x), + Teuchos::rcpFromRef(x_dot), + Teuchos::null, // xdotdot + Teuchos::null, // dxdp + stamp); } void PiroObserver:: -observeSolution(const Thyra_Vector &solution, - const Thyra_MultiVector &solution_dxdp, - const Thyra_Vector &solution_dot, +observeSolution(const Thyra_Vector& x, + const Thyra_MultiVector& dxdp, + const Thyra_Vector& x_dot, const ST stamp) { - this->observeSolutionImpl(solution, solution_dxdp, solution_dot, stamp); - stepper_counter_++; + this->observeSolutionImpl(Teuchos::rcpFromRef(x), + Teuchos::rcpFromRef(x_dot), + Teuchos::null, // xdotdot + Teuchos::rcpFromRef(dxdp), + stamp); } void PiroObserver:: -observeSolution(const Thyra_Vector &solution, - const Thyra_Vector &solution_dot, - const Thyra_Vector &solution_dotdot, +observeSolution(const Thyra_Vector& x, + const Thyra_Vector& x_dot, + const Thyra_Vector& x_dotdot, const ST stamp) { - this->observeSolutionImpl(solution, solution_dot, solution_dotdot, stamp); - stepper_counter_++; + this->observeSolutionImpl(Teuchos::rcpFromRef(x), + Teuchos::rcpFromRef(x_dot), + Teuchos::rcpFromRef(x_dotdot), + Teuchos::null, // dxdp + stamp); } void PiroObserver:: -observeSolution(const Thyra_Vector &solution, - const Thyra_MultiVector &solution_dxdp, - const Thyra_Vector &solution_dot, - const Thyra_Vector &solution_dotdot, +observeSolution(const Thyra_Vector& x, + const Thyra_MultiVector& dxdp, + const Thyra_Vector& x_dot, + const Thyra_Vector& x_dotdot, const ST stamp) { - this->observeSolutionImpl(solution, solution_dxdp, solution_dot, solution_dotdot, stamp); - stepper_counter_++; + this->observeSolutionImpl(Teuchos::rcpFromRef(x), + Teuchos::rcpFromRef(x_dot), + Teuchos::rcpFromRef(x_dotdot), + Teuchos::rcpFromRef(dxdp), + stamp); } void PiroObserver:: -observeSolution(const Thyra_MultiVector &solution, +observeSolution(const Thyra_MultiVector& x, const ST stamp) { - this->observeSolutionImpl(solution, stamp); - stepper_counter_++; -} - -void PiroObserver:: -observeSolution(const Thyra_MultiVector &solution, - const Thyra_MultiVector &solution_dxdp, const ST stamp) -{ - this->observeSolutionImpl(solution, solution_dxdp, stamp); - stepper_counter_++; -} - -void PiroObserver:: -parameterChanged(const std::string& param) -{ - impl_.parameterChanged(param); -} - -void PiroObserver:: -observeSolutionImpl(const Thyra_Vector &solution, - const ST defaultStamp) -{ - // Determine the stamp associated with the snapshot - const ST stamp = impl_.getTimeParamValueOrDefault(defaultStamp); - impl_.observeSolution(stamp, solution, Teuchos::null, - Teuchos::null, Teuchos::null); - - // observe responses - if (observe_responses_ == true) { - if (stepper_counter_ % observe_responses_every_n_steps_ == 0) - this->observeResponse(defaultStamp, Teuchos::rcpFromRef(solution)); - } -} - -void PiroObserver:: -observeSolutionImpl(const Thyra_Vector &solution, - const Thyra_MultiVector &solution_dxdp, - const ST defaultStamp) -{ - // Determine the stamp associated with the snapshot - const ST stamp = impl_.getTimeParamValueOrDefault(defaultStamp); - impl_.observeSolution(stamp, solution, Teuchos::constPtr(solution_dxdp), Teuchos::null, Teuchos::null); - - // observe responses - if (observe_responses_ == true) { - if (stepper_counter_ % observe_responses_every_n_steps_ == 0) - this->observeResponse(defaultStamp, Teuchos::rcpFromRef(solution)); - } -} - -void PiroObserver:: -observeSolutionImpl(const Thyra_Vector &solution, - const Thyra_Vector &solution_dot, - const ST defaultStamp) -{ - // Determine the stamp associated with the snapshot - const ST stamp = impl_.getTimeParamValueOrDefault(defaultStamp); - impl_.observeSolution(stamp, solution, Teuchos::null, - Teuchos::constPtr(solution_dot), Teuchos::null); - - // observe responses - if (observe_responses_ == true) { - if (stepper_counter_ % observe_responses_every_n_steps_ == 0) - this->observeResponse(defaultStamp, - Teuchos::rcpFromRef(solution), - Teuchos::rcpFromRef(solution_dot)); - } + int x_ncols = x.domain()->dim(); + if (x_ncols==1) { + observeSolutionImpl(x.col(0),Teuchos::null,Teuchos::null,Teuchos::null,stamp); + } else if (x_ncols==2) { + observeSolutionImpl(x.col(0),x.col(1),Teuchos::null,Teuchos::null,stamp); + } else { + observeSolutionImpl(x.col(0),x.col(1),x.col(2),Teuchos::null,stamp); + } } void PiroObserver:: -observeSolutionImpl(const Thyra_Vector &solution, - const Thyra_MultiVector &solution_dxdp, - const Thyra_Vector &solution_dot, - const ST defaultStamp) +observeSolution(const Thyra_MultiVector& x, + const Thyra_MultiVector& dxdp, + const ST stamp) { - // Determine the stamp associated with the snapshot - const ST stamp = impl_.getTimeParamValueOrDefault(defaultStamp); - impl_.observeSolution(stamp, solution, Teuchos::constPtr(solution_dxdp), Teuchos::constPtr(solution_dot), Teuchos::null); - - // observe responses - if (observe_responses_ == true) { - if (stepper_counter_ % observe_responses_every_n_steps_ == 0) - this->observeResponse(defaultStamp, - Teuchos::rcpFromRef(solution), - Teuchos::rcpFromRef(solution_dot)); - } + int x_ncols = x.domain()->dim(); + if (x_ncols==1) { + observeSolutionImpl(x.col(0),Teuchos::null,Teuchos::null,Teuchos::rcpFromRef(dxdp),stamp); + } else if (x_ncols==2) { + observeSolutionImpl(x.col(0),x.col(1),Teuchos::null,Teuchos::rcpFromRef(dxdp),stamp); + } else { + observeSolutionImpl(x.col(0),x.col(1),x.col(2),Teuchos::rcpFromRef(dxdp),stamp); + } } void PiroObserver:: -observeSolutionImpl(const Thyra_Vector &solution, - const Thyra_Vector &solution_dot, - const Thyra_Vector &solution_dotdot, +observeSolutionImpl(const Teuchos::RCP& x, + const Teuchos::RCP& x_dot, + const Teuchos::RCP& x_dotdot, + const Teuchos::RCP& dxdp, const ST defaultStamp) { - // Determine the stamp associated with the snapshot - const ST stamp = impl_.getTimeParamValueOrDefault(defaultStamp); - impl_.observeSolution(stamp, solution, Teuchos::null, Teuchos::constPtr(solution_dot), Teuchos::constPtr(solution_dotdot)); - - // observe responses - if (observe_responses_ == true) { - if (stepper_counter_ % observe_responses_every_n_steps_ == 0) - this->observeResponse(defaultStamp, - Teuchos::rcpFromRef(solution), - Teuchos::rcpFromRef(solution_dot), - Teuchos::rcpFromRef(solution_dotdot)); - } -} - + stepper_counter_++; -void PiroObserver:: -observeSolutionImpl(const Thyra_Vector &solution, - const Thyra_MultiVector &solution_dxdp, - const Thyra_Vector &solution_dot, - const Thyra_Vector &solution_dotdot, - const ST defaultStamp) -{ // Determine the stamp associated with the snapshot const ST stamp = impl_.getTimeParamValueOrDefault(defaultStamp); - impl_.observeSolution(stamp, solution, Teuchos::constPtr(solution_dxdp), - Teuchos::constPtr(solution_dot), Teuchos::constPtr(solution_dotdot)); + impl_.observeSolution(stamp, *x, dxdp.ptr(), x_dot.ptr(), x_dotdot.ptr()); // observe responses if (observe_responses_ == true) { if (stepper_counter_ % observe_responses_every_n_steps_ == 0) - this->observeResponse(defaultStamp, - Teuchos::rcpFromRef(solution), - Teuchos::rcpFromRef(solution_dot), - Teuchos::rcpFromRef(solution_dotdot)); + this->observeResponse(defaultStamp,x,x_dot,x_dotdot); } -} - -void PiroObserver:: -observeSolutionImpl(const Thyra_MultiVector &solution, - const ST defaultStamp) -{ - impl_.observeSolution(defaultStamp, solution, Teuchos::null); } -void PiroObserver:: -observeSolutionImpl(const Thyra_MultiVector &solution, - const Thyra_MultiVector &solution_dxdp, - const ST defaultStamp) -{ - impl_.observeSolution(defaultStamp, solution, - Teuchos::constPtr(solution_dxdp)); -} - - void PiroObserver:: observeResponse(const ST defaultStamp, - Teuchos::RCP solution, - Teuchos::RCP solution_dot, - Teuchos::RCP /* solution_dotdot */) + Teuchos::RCP x, + Teuchos::RCP x_dot, + Teuchos::RCP /* x_dotdot */) { - //IKT 5/10/17: note that this function takes solution_dotdot as an input + //IKT 5/10/17: note that this function takes x_dotdot as an input //argument but does not do anything with it yet. This can be modified //if desired. - std::map m_response_index_to_name; - // build out args and evaluate responses if they exist - Thyra::ModelEvaluatorBase::OutArgs outArgs = model_->createOutArgs(); - if(outArgs.Ng()>0) { - // build the in arguments - Thyra::ModelEvaluatorBase::InArgs nominal_values = model_->getNominalValues(); - Thyra::ModelEvaluatorBase::InArgs inArgs = model_->createInArgs(); - inArgs.setArgs(nominal_values); - inArgs.set_x(solution); - if(inArgs.supports(Thyra::ModelEvaluatorBase::IN_ARG_x_dot)) - inArgs.set_x_dot(solution_dot); - if (inArgs.supports(Thyra::ModelEvaluatorBase::IN_ARG_t)) { - const ST time = impl_.getTimeParamValueOrDefault(defaultStamp); - inArgs.set_t(time); - *out << "Time = " << time << "\n"; - } - - // set up the output arguments, in this case only the responses - for(int i=0;iget_g_space(i))); - - // Solve the model - model_->evalModel(inArgs, outArgs); - - std::size_t precision = 8; - std::size_t value_width = precision + 7; - *out << std::scientific << std::showpoint << std::setprecision(precision) << std::left; - + auto outArgs = model_->createOutArgs(); + if (outArgs.Ng()==0) { + return; + } + // build the in arguments + auto nominal_values = model_->getNominalValues(); + auto inArgs = model_->createInArgs(); + inArgs.setArgs(nominal_values); + inArgs.set_x(x); + if (inArgs.supports(Thyra::ModelEvaluatorBase::IN_ARG_x_dot)) + inArgs.set_x_dot(x_dot); + if (inArgs.supports(Thyra::ModelEvaluatorBase::IN_ARG_t)) { + const ST time = impl_.getTimeParamValueOrDefault(defaultStamp); + inArgs.set_t(time); + *out << "Time = " << time << "\n"; + } - // Note that we don't have g_names support in thyra yet. Once - // this is added, we can print response names as well. + // set up the output arguments, in this case only the responses + for(int i=0;iget_g_space(i))); + + // Solve the model + model_->evalModel(inArgs, outArgs); + + std::size_t precision = 8; + std::size_t value_width = precision + 7; + *out << std::scientific << std::showpoint << std::setprecision(precision) << std::left; + + // Note that we don't have g_names support in thyra yet. Once + // this is added, we can print response names as well. + + //OG It seems that outArgs.Ng() always returns 1, so, there is 1 response vector only, Response[0]. + //This response vector contains different responses (min, max, norms) and it would be good + //to have functionality to obtain relative responses only for some values. But it would require more + //parameters in param list. Alternatively, one can rewrite the code below to use is_relative + //as an array of markers for relative responses for Response[0] only. This is not the case + //right now and if in the param list "Relative Responses"="{0}", the code below will compute + //relative values for all terms in vector Response[0]. + if((!firstResponseObtained) && calculateRelativeResponses ){ + storedResponses.resize(outArgs.Ng()); + is_relative.resize(outArgs.Ng(), false); + } - //OG It seems that outArgs.Ng() always returns 1, so, there is 1 response vector only, Response[0]. - //This response vector contains different responses (min, max, norms) and it would be good - //to have functionality to obtain relative responses only for some values. But it would require more - //parameters in param list. Alternatively, one can rewrite the code below to use is_relative - //as an array of markers for relative responses for Response[0] only. This is not the case - //right now and if in the param list "Relative Responses"="{0}", the code below will compute - //relative values for all terms in vector Response[0]. - if((!firstResponseObtained) && calculateRelativeResponses ){ - storedResponses.resize(outArgs.Ng()); - is_relative.resize(outArgs.Ng(), false); + for (int i=0; i::const_iterator itr = m_response_index_to_name.find(i); - if(itr!=m_response_index_to_name.end()) - ss << " Response \"" << itr->second << "\" = "; - else - ss << " Response[" << i << "] = "; - - //ss << "relative resp size? " << relative_responses.size() << "\n"; - //ss << "relative resp values? " << relative_responses << "\n"; - - Teuchos::RCP > g = outArgs.get_g(i); - *out << ss.str(); // " Response[" << i << "] = "; - for(Thyra::Ordinal k=0;kspace()->dim();k++) - *out << std::setw(value_width) << Thyra::get_ele(*g,k) << " "; - *out << std::endl; - - if(firstResponseObtained && calculateRelativeResponses ) - if(is_relative[i]){ - *out << "\n"; - *out << "Relative Response[" << i << "] = "; - for( size_t j = 0; j < storedResponses[i].size(); j++){ - double prevresp = storedResponses[i][j]; - if( std::abs(prevresp) > tol ){ - *out << std::setw(value_width) << (Thyra::get_ele(*g,j) - prevresp)/prevresp << " "; - }else{ - *out << " N/A(int. value 0) "; - } - } - *out << "\n"; + if (not firstResponseObtained and calculateRelativeResponses) { + for(int j = 0; j < relative_responses.size(); j++){ + int resp_index = relative_responses[j]; + if( (resp_index < outArgs.Ng()) ) + is_relative[resp_index] = true; } + } - if( (!firstResponseObtained) && calculateRelativeResponses ){ - for(int j = 0; j < relative_responses.size(); j++){ - int resp_index = relative_responses[j]; - if( (resp_index < outArgs.Ng()) ) - is_relative[resp_index] = true; - } - - } - //Save first responses for relative changes in st - if( (!firstResponseObtained) && calculateRelativeResponses ){ - int gsize = g->space()->dim(); - storedResponses[i].resize(gsize); - for (int j = 0; j < gsize; j++) - storedResponses[i][j] = Thyra::get_ele(*g,j); - }//end if !firstRessponseObtained - }//end of loop over outArgs.Ng() - firstResponseObtained = true; + //Save first responses for relative changes in st + if (not firstResponseObtained and calculateRelativeResponses) { + int gsize = g->space()->dim(); + storedResponses[i].resize(gsize); + for (int j = 0; j < gsize; j++) + storedResponses[i][j] = Thyra::get_ele(*g,j); + } } + firstResponseObtained = true; } } // namespace Albany diff --git a/src/Albany_PiroObserver.hpp b/src/Albany_PiroObserver.hpp index b932450c9a..15f01adc71 100644 --- a/src/Albany_PiroObserver.hpp +++ b/src/Albany_PiroObserver.hpp @@ -24,115 +24,68 @@ class PiroObserver : public Piro::ObserverBase { explicit PiroObserver(const Teuchos::RCP &app, Teuchos::RCP model = Teuchos::null); - virtual void observeSolution( - const Thyra_Vector& solution); - - virtual void observeSolution( - const Thyra_Vector& solution, - const Thyra_MultiVector& solution_dxdp); - - virtual void observeSolution( - const Thyra_Vector& solution, - const ST stamp); + void observeSolution (const Thyra_Vector& x) + { + observeSolution(x,zero()); + } + + void observeSolution (const Thyra_Vector& x, + const Thyra_MultiVector& dxdp) + { + observeSolution(x,dxdp); + } + + void observeSolution (const Thyra_Vector& x, + const ST stamp); + + void observeSolution (const Thyra_Vector& x, + const Thyra_MultiVector& dxdp, + const ST stamp); + + void observeSolution (const Thyra_Vector& x, + const Thyra_Vector& x_dot, + const ST stamp); - virtual void observeSolution( - const Thyra_Vector& solution, - const Thyra_MultiVector& solution_dxdp, - const ST stamp); - - virtual void observeSolution( - const Thyra_Vector& solution, - const Thyra_Vector& solution_dot, - const ST stamp); + void observeSolution (const Thyra_Vector& x, + const Thyra_MultiVector& dxdp, + const Thyra_Vector& x_dot, + const ST stamp); - virtual void observeSolution( - const Thyra_Vector& solution, - const Thyra_MultiVector& solution_dxdp, - const Thyra_Vector& solution_dot, - const ST stamp); + void observeSolution (const Thyra_Vector& x, + const Thyra_Vector& x_dot, + const Thyra_Vector& x_dotdot, + const ST stamp); + + void observeSolution (const Thyra_Vector& x, + const Thyra_MultiVector& dxdp, + const Thyra_Vector& x_dot, + const Thyra_Vector& x_dotdot, + const ST stamp); - virtual void observeSolution( - const Thyra_Vector& solution, - const Thyra_Vector& solution_dot, - const Thyra_Vector& solution_dotdot, - const ST stamp); - - virtual void observeSolution( - const Thyra_Vector& solution, - const Thyra_MultiVector& solution_dxdp, - const Thyra_Vector& solution_dot, - const Thyra_Vector& solution_dotdot, - const ST stamp); + void observeSolution (const Thyra_MultiVector& solution, + const ST stamp); - virtual void observeSolution( - const Thyra_MultiVector& solution, - const ST stamp); - - virtual void observeSolution( - const Thyra_MultiVector& solution, - const Thyra_MultiVector& solution_dxdp, - const ST stamp); + void observeSolution (const Thyra_MultiVector& solution, + const Thyra_MultiVector& solution_dxdp, + const ST stamp); - virtual void parameterChanged( - const std::string& param); + void parameterChanged (const std::string& param) { impl_.parameterChanged(param); } private: - void observeSolutionImpl( - const Thyra_Vector& solution, - const ST defaultStamp); - - void observeSolutionImpl( - const Thyra_Vector& solution, - const Thyra_MultiVector& solution_dxdp, - const ST defaultStamp); - - void observeSolutionImpl( - const Thyra_Vector& solution, - const Thyra_Vector& solution_dot, - const ST defaultStamp); - - void observeSolutionImpl( - const Thyra_Vector& solution, - const Thyra_MultiVector& solution_dxdp, - const Thyra_Vector& solution_dot, - const ST defaultStamp); - - void observeSolutionImpl( - const Thyra_Vector& solution, - const Thyra_Vector& solution_dot, - const Thyra_Vector& solution_dotdot, - const ST defaultStamp); - - void observeSolutionImpl( - const Thyra_Vector& solution, - const Thyra_MultiVector& solution_dxdp, - const Thyra_Vector& solution_dot, - const Thyra_Vector& solution_dotdot, - const ST defaultStamp); - - void observeSolutionImpl( - const Thyra_MultiVector &solution, - const ST defaultStamp); - - void observeSolutionImpl( - const Thyra_MultiVector &solution, - const Thyra_MultiVector &solution_dxdp, - const ST defaultStamp); - - void observeTpetraSolutionImpl( - const Tpetra_Vector &solution, - Teuchos::Ptr solution_dot, - Teuchos::Ptr solution_dotdot, - const ST defaultStamp); + + void observeSolutionImpl (const Teuchos::RCP& x, + const Teuchos::RCP& x_dot, + const Teuchos::RCP& x_dotdot, + const Teuchos::RCP& dxdp, + const ST defaultStamp); // The following function is for calculating / printing responses every step. // It is currently not implemented for the case of an Teuchos::RCP // argument; this may be desired at some point in the future. - void observeResponse( - const ST defaultStamp, - Teuchos::RCP solution, - Teuchos::RCP solution_dot = Teuchos::null, - Teuchos::RCP solution_dotdot = Teuchos::null); + void observeResponse(const ST defaultStamp, + Teuchos::RCP x, + Teuchos::RCP x_dot, + Teuchos::RCP x_dotdot); ObserverImpl impl_; @@ -140,6 +93,8 @@ class PiroObserver : public Piro::ObserverBase { protected: + static ST zero () { return Teuchos::ScalarTraits::zero(); } + bool observe_responses_; int stepper_counter_; diff --git a/src/Albany_SolverFactory.cpp b/src/Albany_SolverFactory.cpp index 5e963e2264..28dc854e20 100644 --- a/src/Albany_SolverFactory.cpp +++ b/src/Albany_SolverFactory.cpp @@ -155,9 +155,8 @@ SolverFactory::createModel (const Teuchos::RCP& app, Teuchos::RCP> SolverFactory:: -createSolver (const Teuchos::RCP& solverComm, - const Teuchos::RCP& model_tmp, - const Teuchos::RCP& adjointModel_tmp, +createSolver (const Teuchos::RCP& model_tmp, + const Teuchos::RCP& adjointModel_tmp, const bool forwardMode) { const auto piroParams = Teuchos::sublist(m_appParams, "Piro"); @@ -254,24 +253,13 @@ createSolver (const Teuchos::RCP& solverComm, } } - const auto app = model_tmp->getAlbanyApp(); - const auto solMgr = app->getAdaptSolMgr(); + const auto app = model_tmp->getAlbanyApp(); - Piro::SolverFactory piroFactory; - m_observer = Teuchos::rcp(new PiroObserver(app, modelWithSolve)); + auto observer = Teuchos::rcp(new PiroObserver(app, modelWithSolve)); + Piro::SolverFactory piroFactory; return piroFactory.createSolver( - piroParams, modelWithSolve, adjointModelWithSolve, m_observer); - TEUCHOS_TEST_FOR_EXCEPTION( - true, - std::logic_error, - "Reached end of createModel()" - << "\n"); - - // Silence compiler warning in case it wasn't used (due to ifdef logic) - (void) solverComm; - - return Teuchos::null; + piroParams, modelWithSolve, adjointModelWithSolve, observer); } void SolverFactory:: diff --git a/src/Albany_SolverFactory.hpp b/src/Albany_SolverFactory.hpp index 63f9d6456b..b4289fb285 100644 --- a/src/Albany_SolverFactory.hpp +++ b/src/Albany_SolverFactory.hpp @@ -50,9 +50,8 @@ class SolverFactory { const bool adjoint_model = false); Teuchos::RCP> - createSolver (const Teuchos::RCP& solverComm, - const Teuchos::RCP& model, - const Teuchos::RCP& adjointModel, + createSolver (const Teuchos::RCP& model, + const Teuchos::RCP& adjointModel, const bool forwardMode); Teuchos::ParameterList& @@ -67,12 +66,6 @@ class SolverFactory { return m_appParams; } - Teuchos::RCP> - returnObserver() const - { - return m_observer; - }; - // Functions to generate reference parameter lists for validation // EGN 9/2013: made these three functions public, as they pertain to valid // parameter lists for Albany::Application objects, which may get created @@ -92,8 +85,6 @@ class SolverFactory { void setSolverParamDefaults(Teuchos::ParameterList* appParams, int myRank); - Teuchos::RCP> m_observer; - //! Parameter list specifying what solver to create Teuchos::RCP m_appParams; diff --git a/src/Albany_StatelessObserverImpl.cpp b/src/Albany_StatelessObserverImpl.cpp index 102ab1b768..8cfc406e16 100644 --- a/src/Albany_StatelessObserverImpl.cpp +++ b/src/Albany_StatelessObserverImpl.cpp @@ -36,98 +36,33 @@ getTimeParamValueOrDefault (RealType defaultValue) const { return this_time; } -Teuchos::RCP -StatelessObserverImpl::getNonOverlappedVectorSpace () const { - return app_->getVectorSpace(); -} - -void StatelessObserverImpl::observeSolution ( - double stamp, - const Thyra_Vector &nonOverlappedSolution, - const Teuchos::Ptr &nonOverlappedSolution_dxdp, - const Teuchos::Ptr& nonOverlappedSolutionDot) -{ - Teuchos::TimeMonitor timer(*solOutTime_); - const Teuchos::RCP overlappedSolution = - app_->getAdaptSolMgr()->updateAndReturnOverlapSolution(nonOverlappedSolution); - Teuchos::RCP overlappedSolutionDxDp = Teuchos::null; - if (nonOverlappedSolution_dxdp != Teuchos::null) { - overlappedSolutionDxDp = app_->getAdaptSolMgr()->updateAndReturnOverlapSolutionDxDp(*nonOverlappedSolution_dxdp); - /*auto out_ = Teuchos::VerboseObjectBase::getDefaultOStream(); - int num_param = overlappedSolutionDxDp->domain()->dim(); - for (int np = 0; np < num_param; np++) { - *out_ << "\n*** StatelessObserverImpl::observeSolution overlappedSolutionDxDp" << np << " ***\n"; - Teuchos::RCP> solution_dxdp_np = overlappedSolutionDxDp->col(np); - Teuchos::Range1D range; - RTOpPack::ConstSubVectorView dxdpv; - solution_dxdp_np->acquireDetachedView(range, &dxdpv); - auto dxdpa = dxdpv.values(); - for (auto i = 0; i < dxdpa.size(); ++i) *out_ << dxdpa[i] << " "; - *out_ << "\n*** StatelessObserverImpl::observeSolution overlappedSolutionDxDp" << np << " ***\n"; - }*/ - } - if (nonOverlappedSolutionDot != Teuchos::null) { - const Teuchos::RCP overlappedSolutionDot = - app_->getAdaptSolMgr()->updateAndReturnOverlapSolutionDot(*nonOverlappedSolutionDot); - app_->getDiscretization()->writeSolution( - *overlappedSolution, overlappedSolutionDxDp, *overlappedSolutionDot, stamp, true, - force_write_solution_); - } - else { - app_->getDiscretization()->writeSolution(*overlappedSolution, overlappedSolutionDxDp, - stamp, true, force_write_solution_); - } -} - -void StatelessObserverImpl::observeSolution ( - double stamp, - const Thyra_Vector &nonOverlappedSolution, - const Teuchos::Ptr &nonOverlappedSolution_dxdp, - const Teuchos::Ptr& nonOverlappedSolutionDot, - const Teuchos::Ptr& nonOverlappedSolutionDotDot) +void StatelessObserverImpl:: +observeSolution (double stamp, + const Teuchos::RCP& x, + const Teuchos::RCP& x_dot, + const Teuchos::RCP& x_dotdot, + const Teuchos::RCP& dxdp) { Teuchos::TimeMonitor timer(*solOutTime_); - const Teuchos::RCP overlappedSolution = - app_->getAdaptSolMgr()->updateAndReturnOverlapSolution(nonOverlappedSolution); - Teuchos::RCP overlappedSolutionDxDp = Teuchos::null; - if (nonOverlappedSolution_dxdp != Teuchos::null) { - overlappedSolutionDxDp = app_->getAdaptSolMgr()->updateAndReturnOverlapSolutionDxDp(*nonOverlappedSolution_dxdp); + auto overlapped_x = app_->getAdaptSolMgr()->updateAndReturnOverlapSolution(*x); + Teuchos::RCP overlapped_dxdp; + if (dxdp != Teuchos::null) { + overlapped_dxdp = app_->getAdaptSolMgr()->updateAndReturnOverlapSolutionDxDp(*dxdp); } - if (nonOverlappedSolutionDot != Teuchos::null) { - const Teuchos::RCP overlappedSolutionDot = - app_->getAdaptSolMgr()->updateAndReturnOverlapSolutionDot(*nonOverlappedSolutionDot); - if (nonOverlappedSolutionDotDot != Teuchos::null) { - const Teuchos::RCP overlappedSolutionDotDot = - app_->getAdaptSolMgr()->updateAndReturnOverlapSolutionDotDot(*nonOverlappedSolutionDotDot); + if (x_dot != Teuchos::null) { + auto overlapped_x_dot = app_->getAdaptSolMgr()->updateAndReturnOverlapSolutionDot(*x_dot); + if (x_dotdot != Teuchos::null) { + auto overlapped_x_dotdot = app_->getAdaptSolMgr()->updateAndReturnOverlapSolutionDotDot(*x_dotdot); app_->getDiscretization()->writeSolution( - *overlappedSolution, overlappedSolutionDxDp, *overlappedSolutionDot, *overlappedSolutionDotDot, + *overlapped_x, overlapped_dxdp, *overlapped_x_dot, *overlapped_x_dotdot, stamp, true, force_write_solution_); - } - else { + } else { app_->getDiscretization()->writeSolution( - *overlappedSolution, overlappedSolutionDxDp, *overlappedSolutionDot, - stamp, true, force_write_solution_); + *overlapped_x, overlapped_dxdp, *overlapped_x_dot, stamp, true, force_write_solution_); } + } else { + app_->getDiscretization()->writeSolution(*overlapped_x, overlapped_dxdp, stamp, true, force_write_solution_); } - else { - app_->getDiscretization()->writeSolution( - *overlappedSolution, overlappedSolutionDxDp, stamp, true, force_write_solution_); - } -} - -void StatelessObserverImpl::observeSolution ( - double stamp, const Thyra_MultiVector &nonOverlappedSolution, - const Teuchos::Ptr &nonOverlappedSolution_dxdp) -{ - Teuchos::TimeMonitor timer(*solOutTime_); - const Teuchos::RCP overlappedSolution = - app_->getAdaptSolMgr()->updateAndReturnOverlapSolutionMV(nonOverlappedSolution); - Teuchos::RCP overlappedSolutionDxDp = Teuchos::null; - if (nonOverlappedSolution_dxdp != Teuchos::null) { - overlappedSolutionDxDp = app_->getAdaptSolMgr()->updateAndReturnOverlapSolutionDxDp(*nonOverlappedSolution_dxdp); - } - app_->getDiscretization()->writeSolutionMV(*overlappedSolution, overlappedSolutionDxDp, - stamp, true, force_write_solution_); } } // namespace Albany diff --git a/src/Albany_StatelessObserverImpl.hpp b/src/Albany_StatelessObserverImpl.hpp index b8615fca7c..536f2cc963 100644 --- a/src/Albany_StatelessObserverImpl.hpp +++ b/src/Albany_StatelessObserverImpl.hpp @@ -27,25 +27,11 @@ class StatelessObserverImpl { RealType getTimeParamValueOrDefault(RealType defaultValue) const; - Teuchos::RCP getNonOverlappedVectorSpace() const; - - virtual void observeSolution ( - double stamp, - const Thyra_Vector& nonOverlappedSolution, - const Teuchos::Ptr& nonOverlappedSolution_dxdp, - const Teuchos::Ptr& nonOverlappedSolutionDot, - const Teuchos::Ptr& nonOverlappedSolutionDotDot); - - virtual void observeSolution ( - double stamp, - const Thyra_Vector& nonOverlappedSolution, - const Teuchos::Ptr& nonOverlappedSolution_dxdp, - const Teuchos::Ptr& nonOverlappedSolutionDot); - - virtual void observeSolution ( - double stamp, - const Thyra_MultiVector& nonOverlappedSolution, - const Teuchos::Ptr& nonOverlappedSolution_dxdp); + void observeSolution (double stamp, + const Teuchos::RCP& x, + const Teuchos::RCP& x_dot, + const Teuchos::RCP& x_dotdot, + const Teuchos::RCP& x_dxdp); protected: Teuchos::RCP app_; diff --git a/src/Main_Analysis.cpp b/src/Main_Analysis.cpp index 1651ced4f9..f008b361b5 100644 --- a/src/Main_Analysis.cpp +++ b/src/Main_Analysis.cpp @@ -81,7 +81,7 @@ int main(int argc, char *argv[]) { && slvrfctry.getParameters()->sublist("Piro").sublist("Analysis").get("Transient"); const auto albanyAdjointModel = explicitMatrixTranspose || transientAnalysis ? slvrfctry.createModel(albanyApp, true) : Teuchos::null; - const auto solver = slvrfctry.createSolver(comm, albanyModel, albanyAdjointModel, false); + const auto solver = slvrfctry.createSolver(albanyModel, albanyAdjointModel, false); stackedTimer->stop("Albany: Setup Time"); diff --git a/src/Main_Solve.cpp b/src/Main_Solve.cpp index 7d9faa942e..aa44e8e4f4 100644 --- a/src/Main_Solve.cpp +++ b/src/Main_Solve.cpp @@ -139,7 +139,7 @@ int main(int argc, char *argv[]) // Explicit adjoint model is not needed if we are not computing adjoint sensitivities const bool explicitAdjointModel = albanyApp->isAdjointSensitivities() && explicitMatrixTranspose; const auto albanyAdjointModel = explicitAdjointModel ? slvrfctry.createModel(albanyApp, true) : Teuchos::null; - const auto solver = slvrfctry.createSolver(comm, albanyModel, albanyAdjointModel, true); + const auto solver = slvrfctry.createSolver(albanyModel, albanyAdjointModel, true); stackedTimer->stop("Albany: Setup Time"); diff --git a/src/landIce/interfaceWithMPAS/Interface.cpp b/src/landIce/interfaceWithMPAS/Interface.cpp index cff5499ae7..d5cfcd1f2b 100644 --- a/src/landIce/interfaceWithMPAS/Interface.cpp +++ b/src/landIce/interfaceWithMPAS/Interface.cpp @@ -315,7 +315,7 @@ void velocity_solver_solve_fo(int nLayers, int globalVerticesStride, Teuchos::ArrayRCP solution_constView; try { auto model = slvrfctry->createModel(albanyApp); - solver = slvrfctry->createSolver(mpiComm, model, Teuchos::null, true); + solver = slvrfctry->createSolver(model, Teuchos::null, true); Teuchos::ParameterList solveParams; solveParams.set("Compute Sensitivities", false);