Panzer Version of the Day
Loading...
Searching...
No Matches
Panzer_ModelEvaluator_impl.hpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Panzer: A partial differential equation assembly
5// engine for strongly coupled complex multiphysics systems
6// Copyright (2011) Sandia Corporation
7//
8// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9// the U.S. Government retains certain rights in this software.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact Roger P. Pawlowski (rppawlo@sandia.gov) and
39// Eric C. Cyr (eccyr@sandia.gov)
40// ***********************************************************************
41// @HEADER
42
43#ifndef __Panzer_ModelEvaluator_impl_hpp__
44#define __Panzer_ModelEvaluator_impl_hpp__
45
46#include "Teuchos_DefaultComm.hpp"
47#include "Teuchos_ArrayRCP.hpp"
48
49#include "PanzerDiscFE_config.hpp"
50#include "Panzer_Traits.hpp"
57#include "Panzer_GlobalData.hpp"
65
66#include "Thyra_TpetraThyraWrappers.hpp"
67#include "Thyra_SpmdVectorBase.hpp"
68#include "Thyra_DefaultSpmdVector.hpp"
69#include "Thyra_DefaultSpmdVectorSpace.hpp"
70#include "Thyra_DefaultMultiVectorProductVectorSpace.hpp"
71#include "Thyra_DefaultMultiVectorProductVector.hpp"
72#include "Thyra_MultiVectorStdOps.hpp"
73#include "Thyra_VectorStdOps.hpp"
74
75// For writing out residuals/Jacobians
76#include "Thyra_ProductVectorBase.hpp"
77#include "Thyra_BlockedLinearOpBase.hpp"
78#include "Thyra_TpetraVector.hpp"
79#include "Thyra_TpetraLinearOp.hpp"
80#include "Tpetra_CrsMatrix.hpp"
81
82// Constructors/Initializers/Accessors
83
84template<typename Scalar>
86ModelEvaluator(const Teuchos::RCP<panzer::FieldManagerBuilder>& fmb,
87 const Teuchos::RCP<panzer::ResponseLibrary<panzer::Traits> >& rLibrary,
88 const Teuchos::RCP<const panzer::LinearObjFactory<panzer::Traits> >& lof,
89 const std::vector<Teuchos::RCP<Teuchos::Array<std::string> > >& p_names,
90 const std::vector<Teuchos::RCP<Teuchos::Array<double> > >& p_values,
91 const Teuchos::RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> > & solverFactory,
92 const Teuchos::RCP<panzer::GlobalData>& global_data,
93 bool build_transient_support,
94 double t_init)
95 : t_init_(t_init)
96 , num_me_parameters_(0)
97 , do_fd_dfdp_(false)
98 , fd_perturb_size_(1e-7)
99 , require_in_args_refresh_(true)
100 , require_out_args_refresh_(true)
101 , responseLibrary_(rLibrary)
102 , global_data_(global_data)
103 , build_transient_support_(build_transient_support)
104 , lof_(lof)
105 , solverFactory_(solverFactory)
106 , oneTimeDirichletBeta_on_(false)
107 , oneTimeDirichletBeta_(0.0)
108 , build_volume_field_managers_(true)
109 , build_bc_field_managers_(true)
110 , active_evaluation_types_(Sacado::mpl::size<panzer::Traits::EvalTypes>::value, true)
111 , write_matrix_count_(0)
112{
113 using Teuchos::RCP;
114 using Teuchos::rcp;
115 using Teuchos::rcp_dynamic_cast;
116 using Teuchos::tuple;
117 using Thyra::VectorBase;
118 using Thyra::createMember;
119
120 TEUCHOS_ASSERT(lof_!=Teuchos::null);
121
123 ae_tm_.buildObjects(builder);
124
125 //
126 // Build x, f spaces
127 //
128
129 // dynamic cast to blocked LOF for now
130 RCP<const ThyraObjFactory<Scalar> > tof = rcp_dynamic_cast<const ThyraObjFactory<Scalar> >(lof,true);
131
132 x_space_ = tof->getThyraDomainSpace();
133 f_space_ = tof->getThyraRangeSpace();
134
135 //
136 // Setup parameters
137 //
138 for(std::size_t i=0;i<p_names.size();i++)
139 addParameter(*(p_names[i]),*(p_values[i]));
140
141 // now that the vector spaces are setup we can allocate the nominal values
142 // (i.e. initial conditions)
144}
145
146template<typename Scalar>
148ModelEvaluator(const Teuchos::RCP<const panzer::LinearObjFactory<panzer::Traits> >& lof,
149 const Teuchos::RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> > & solverFactory,
150 const Teuchos::RCP<panzer::GlobalData>& global_data,
151 bool build_transient_support,double t_init)
152 : t_init_(t_init)
153 , num_me_parameters_(0)
154 , do_fd_dfdp_(false)
155 , fd_perturb_size_(1e-7)
156 , require_in_args_refresh_(true)
157 , require_out_args_refresh_(true)
158 , global_data_(global_data)
159 , build_transient_support_(build_transient_support)
160 , lof_(lof)
161 , solverFactory_(solverFactory)
162 , oneTimeDirichletBeta_on_(false)
163 , oneTimeDirichletBeta_(0.0)
164 , build_volume_field_managers_(true)
165 , build_bc_field_managers_(true)
166 , active_evaluation_types_(Sacado::mpl::size<panzer::Traits::EvalTypes>::value, true)
167 , write_matrix_count_(0)
168{
169 using Teuchos::RCP;
170 using Teuchos::rcp_dynamic_cast;
171
172 TEUCHOS_ASSERT(lof_!=Teuchos::null);
173
174 //
175 // Build x, f spaces
176 //
177
178 // dynamic cast to blocked LOF for now
179 RCP<const ThyraObjFactory<Scalar> > tof = rcp_dynamic_cast<const ThyraObjFactory<Scalar> >(lof_,true);
180
181 x_space_ = tof->getThyraDomainSpace();
182 f_space_ = tof->getThyraRangeSpace();
183
184 // now that the vector spaces are setup we can allocate the nominal values
185 // (i.e. initial conditions)
187
188 // allocate a response library so that responses can be added, it will be initialized in "setupModel"
190}
191
192template<typename Scalar>
195{
196 TEUCHOS_ASSERT(false);
197}
198
199// Public functions overridden from ModelEvaulator
200
201template<typename Scalar>
202Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
204{
205 return x_space_;
206}
207
208
209template<typename Scalar>
210Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
212{
213 return f_space_;
214}
215
216template<typename Scalar>
217Teuchos::RCP<const Teuchos::Array<std::string> >
219{
220 TEUCHOS_TEST_FOR_EXCEPTION(!(i>=0 && i<num_me_parameters_),std::runtime_error,
221 "panzer::ModelEvaluator::get_p_names: Requested parameter index out of range.");
222
223 if (i < Teuchos::as<int>(parameters_.size()))
224 return parameters_[i]->names;
225 else if (i < Teuchos::as<int>(parameters_.size()+tangent_space_.size())) {
226 Teuchos::RCP<Teuchos::Array<std::string> > names = rcp(new Teuchos::Array<std::string>);
227 int param_index = i-parameters_.size();
228 std::ostringstream ss;
229 ss << "TANGENT VECTOR: " << param_index;
230 names->push_back(ss.str());
231 return names;
232 }
233 else if (build_transient_support_ && i < Teuchos::as<int>(parameters_.size()+2*tangent_space_.size())) {
234 Teuchos::RCP<Teuchos::Array<std::string> > names = rcp(new Teuchos::Array<std::string>);
235 int param_index = i-parameters_.size()-tangent_space_.size();
236 std::ostringstream ss;
237 ss << "TIME DERIVATIVE TANGENT VECTOR: " << param_index;
238 names->push_back(ss.str());
239 return names;
240 }
241
242 return Teuchos::null;
243}
244
245template<typename Scalar>
246Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
248{
249 TEUCHOS_TEST_FOR_EXCEPTION(!(i>=0 && i<num_me_parameters_),std::runtime_error,
250 "panzer::ModelEvaluator::get_p_space: Requested parameter index out of range.");
251
252 if (i < Teuchos::as<int>(parameters_.size()))
253 return parameters_[i]->space;
254 else if (i < Teuchos::as<int>(parameters_.size()+tangent_space_.size()))
255 return tangent_space_[i-parameters_.size()];
256 else if (build_transient_support_ && i < Teuchos::as<int>(parameters_.size()+2*tangent_space_.size()))
257 return tangent_space_[i-parameters_.size()-tangent_space_.size()];
258
259 return Teuchos::null;
260}
261
262template<typename Scalar>
263Teuchos::ArrayView<const std::string>
265{
266 TEUCHOS_TEST_FOR_EXCEPTION(!(i>=0 && i<Teuchos::as<int>(responses_.size())),std::runtime_error,
267 "panzer::ModelEvaluator::get_g_names: Requested response index out of range.");
268
269 return Teuchos::ArrayView<const std::string>(&(responses_[i]->name),1);
270}
271
272template<typename Scalar>
273const std::string &
275{
276 TEUCHOS_ASSERT(i>=0 &&
277 static_cast<typename std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > >::size_type>(i)<responses_.size());
278
279 return responses_[i]->name;
280}
281
282template<typename Scalar>
283Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
285{
286 TEUCHOS_ASSERT(i>=0 &&
287 static_cast<typename std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > >::size_type>(i)<responses_.size());
288
289 return responses_[i]->space;
290}
291
292template<typename Scalar>
293Thyra::ModelEvaluatorBase::InArgs<Scalar>
295{
296 return getNominalValues();
297}
298
299template<typename Scalar>
300Thyra::ModelEvaluatorBase::InArgs<Scalar>
302{
303 using Teuchos::RCP;
304 using Teuchos::rcp_dynamic_cast;
305
306 if(require_in_args_refresh_) {
307 typedef Thyra::ModelEvaluatorBase MEB;
308
309 //
310 // Refresh nominal values, this is primarily adding
311 // new parameters.
312 //
313
314 MEB::InArgsSetup<Scalar> nomInArgs;
315 nomInArgs = nominalValues_;
316 nomInArgs.setSupports(nominalValues_);
317
318 // setup parameter support
319 nomInArgs.set_Np(num_me_parameters_);
320 for(std::size_t p=0;p<parameters_.size();p++) {
321 // setup nominal in arguments
322 nomInArgs.set_p(p,parameters_[p]->initial_value);
323
324 // We explicitly do not set nominal values for tangent parameters
325 // as these are parameters that should be hidden from client code
326 }
327
328 nominalValues_ = nomInArgs;
329 }
330
331 // refresh no longer required
332 require_in_args_refresh_ = false;
333
334 return nominalValues_;
335}
336
337template<typename Scalar>
338void
340{
341 typedef Thyra::ModelEvaluatorBase MEB;
342
343 //
344 // Setup nominal values
345 //
346
347 MEB::InArgsSetup<Scalar> nomInArgs;
348 nomInArgs.setModelEvalDescription(this->description());
349 nomInArgs.setSupports(MEB::IN_ARG_x);
350 Teuchos::RCP<Thyra::VectorBase<Scalar> > x_nom = Thyra::createMember(x_space_);
351 Thyra::assign(x_nom.ptr(),0.0);
352 nomInArgs.set_x(x_nom);
353 if(build_transient_support_) {
354 nomInArgs.setSupports(MEB::IN_ARG_x_dot,true);
355 nomInArgs.setSupports(MEB::IN_ARG_t,true);
356 nomInArgs.setSupports(MEB::IN_ARG_alpha,true);
357 nomInArgs.setSupports(MEB::IN_ARG_beta,true);
358 nomInArgs.setSupports(MEB::IN_ARG_step_size,true);
359 nomInArgs.setSupports(MEB::IN_ARG_stage_number,true);
360
361 Teuchos::RCP<Thyra::VectorBase<Scalar> > x_dot_nom = Thyra::createMember(x_space_);
362 Thyra::assign(x_dot_nom.ptr(),0.0);
363 nomInArgs.set_x_dot(x_dot_nom);
364 nomInArgs.set_t(t_init_);
365 nomInArgs.set_alpha(0.0); // these have no meaning initially!
366 nomInArgs.set_beta(0.0);
367 //TODO: is this needed?
368 nomInArgs.set_step_size(0.0);
369 nomInArgs.set_stage_number(1.0);
370 }
371
372 // setup parameter support -- for each scalar parameter we support the parameter itself and tangent vectors for x, xdot
373 nomInArgs.set_Np(num_me_parameters_);
374 std::size_t v_index = 0;
375 for(std::size_t p=0;p<parameters_.size();p++) {
376 nomInArgs.set_p(p,parameters_[p]->initial_value);
377 if (!parameters_[p]->is_distributed) {
378 Teuchos::RCP<Thyra::VectorBase<Scalar> > v_nom_x = Thyra::createMember(*tangent_space_[v_index]);
379 Thyra::assign(v_nom_x.ptr(),0.0);
380 nomInArgs.set_p(v_index+parameters_.size(),v_nom_x);
381 if (build_transient_support_) {
382 Teuchos::RCP<Thyra::VectorBase<Scalar> > v_nom_xdot = Thyra::createMember(*tangent_space_[v_index]);
383 Thyra::assign(v_nom_xdot.ptr(),0.0);
384 nomInArgs.set_p(v_index+parameters_.size()+tangent_space_.size(),v_nom_xdot);
385 }
386 ++v_index;
387 }
388 }
389
390 nominalValues_ = nomInArgs;
391}
392
393template <typename Scalar>
395buildVolumeFieldManagers(const bool value)
396{
397 build_volume_field_managers_ = value;
398}
399
400template <typename Scalar>
402buildBCFieldManagers(const bool value)
403{
404 build_bc_field_managers_ = value;
405}
406
407template <typename Scalar>
409setupModel(const Teuchos::RCP<panzer::WorksetContainer> & wc,
410 const std::vector<Teuchos::RCP<panzer::PhysicsBlock> >& physicsBlocks,
411 const std::vector<panzer::BC> & bcs,
412 const panzer::EquationSetFactory & eqset_factory,
413 const panzer::BCStrategyFactory& bc_factory,
416 const Teuchos::ParameterList& closure_models,
417 const Teuchos::ParameterList& user_data,
418 bool writeGraph,const std::string & graphPrefix,
419 const Teuchos::ParameterList& me_params)
420{
421 // First: build residual assembly engine
423 PANZER_FUNC_TIME_MONITOR_DIFF("panzer::ModelEvaluator::setupModel()",setupModel);
424
425 {
426 // 1. build Field manager builder
428
429 Teuchos::RCP<panzer::FieldManagerBuilder> fmb;
430 {
431 PANZER_FUNC_TIME_MONITOR_DIFF("allocate FieldManagerBuilder",allocFMB);
432 fmb = Teuchos::rcp(new panzer::FieldManagerBuilder);
433 fmb->setActiveEvaluationTypes(active_evaluation_types_);
434 }
435 {
436 PANZER_FUNC_TIME_MONITOR_DIFF("fmb->setWorksetContainer()",setupWorksets);
437 fmb->setWorksetContainer(wc);
438 }
439 if (build_volume_field_managers_) {
440 PANZER_FUNC_TIME_MONITOR_DIFF("fmb->setupVolumeFieldManagers()",setupVolumeFieldManagers);
441 fmb->setupVolumeFieldManagers(physicsBlocks,volume_cm_factory,closure_models,*lof_,user_data);
442 }
443 if (build_bc_field_managers_) {
444 PANZER_FUNC_TIME_MONITOR_DIFF("fmb->setupBCFieldManagers()",setupBCFieldManagers);
445 fmb->setupBCFieldManagers(bcs,physicsBlocks,eqset_factory,bc_cm_factory,bc_factory,closure_models,*lof_,user_data);
446 }
447
448 // Print Phalanx DAGs
449 if (writeGraph){
450 if (build_volume_field_managers_)
451 fmb->writeVolumeGraphvizDependencyFiles(graphPrefix, physicsBlocks);
452 if (build_bc_field_managers_)
453 fmb->writeBCGraphvizDependencyFiles(graphPrefix+"BC_");
454 }
455
456 {
457 PANZER_FUNC_TIME_MONITOR_DIFF("AssemblyEngine_TemplateBuilder::buildObjects()",AETM_BuildObjects);
459 ae_tm_.buildObjects(builder);
460 }
461 }
462
463 // Second: build the responses
465
466 {
467 PANZER_FUNC_TIME_MONITOR_DIFF("build response library",buildResponses);
468
469 responseLibrary_->initialize(wc,lof_->getRangeGlobalIndexer(),lof_);
470
471 buildResponses(physicsBlocks,eqset_factory,volume_cm_factory,closure_models,user_data,writeGraph,graphPrefix+"Responses_");
472 buildDistroParamDfDp_RL(wc,physicsBlocks,bcs,eqset_factory,bc_factory,volume_cm_factory,closure_models,user_data,writeGraph,graphPrefix+"Response_DfDp_");
473 buildDistroParamDgDp_RL(wc,physicsBlocks,bcs,eqset_factory,bc_factory,volume_cm_factory,closure_models,user_data,writeGraph,graphPrefix+"Response_DgDp_");
474
475 do_fd_dfdp_ = false;
476 fd_perturb_size_ = 1.0e-7;
477 if (me_params.isParameter("FD Forward Sensitivities"))
478 do_fd_dfdp_ = me_params.get<bool>("FD Forward Sensitivities");
479 if (me_params.isParameter("FD Perturbation Size"))
480 fd_perturb_size_ = me_params.get<double>("FD Perturbation Size");
481 }
482}
483
484template <typename Scalar>
486setupAssemblyInArgs(const Thyra::ModelEvaluatorBase::InArgs<Scalar> & inArgs,
487 panzer::AssemblyEngineInArgs & ae_inargs) const
488{
489 using Teuchos::RCP;
490 using Teuchos::rcp;
491 using Teuchos::rcp_dynamic_cast;
492 using Teuchos::rcp_const_cast;
493 typedef Thyra::ModelEvaluatorBase MEB;
494
495 // if neccessary build a ghosted container
496 if(Teuchos::is_null(ghostedContainer_)) {
497 ghostedContainer_ = lof_->buildGhostedLinearObjContainer();
498 lof_->initializeGhostedContainer(panzer::LinearObjContainer::X |
501 panzer::LinearObjContainer::Mat, *ghostedContainer_);
502 }
503
504 bool is_transient = false;
505 if (inArgs.supports(MEB::IN_ARG_x_dot ))
506 is_transient = !Teuchos::is_null(inArgs.get_x_dot());
507
508 if(Teuchos::is_null(xContainer_))
509 xContainer_ = lof_->buildReadOnlyDomainContainer();
510 if(Teuchos::is_null(xdotContainer_) && is_transient)
511 xdotContainer_ = lof_->buildReadOnlyDomainContainer();
512
513 const RCP<const Thyra::VectorBase<Scalar> > x = inArgs.get_x();
514 RCP<const Thyra::VectorBase<Scalar> > x_dot; // possibly empty, but otherwise uses x_dot
515
516 // Make sure construction built in transient support
517 TEUCHOS_TEST_FOR_EXCEPTION(is_transient && !build_transient_support_, std::runtime_error,
518 "ModelEvaluator was not built with transient support enabled!");
519
520 ae_inargs.container_ = lof_->buildLinearObjContainer(); // we use a new global container
521 ae_inargs.ghostedContainer_ = ghostedContainer_; // we can reuse the ghosted container
522 ae_inargs.alpha = 0.0;
523 ae_inargs.beta = 1.0;
524 ae_inargs.evaluate_transient_terms = false;
525 if (build_transient_support_) {
526 x_dot = inArgs.get_x_dot();
527 ae_inargs.alpha = inArgs.get_alpha();
528 ae_inargs.beta = inArgs.get_beta();
529 ae_inargs.time = inArgs.get_t();
530
531 ae_inargs.step_size= inArgs.get_step_size();
532 ae_inargs.stage_number = inArgs.get_stage_number();
533 ae_inargs.evaluate_transient_terms = true;
534 }
535
536 // this member is handled in the individual functions
537 ae_inargs.apply_dirichlet_beta = false;
538
539 // Set input parameters
540 int num_param_vecs = parameters_.size();
541 for (int i=0; i<num_param_vecs; i++) {
542
543 RCP<const Thyra::VectorBase<Scalar> > paramVec = inArgs.get_p(i);
544 if ( paramVec!=Teuchos::null && !parameters_[i]->is_distributed) {
545 // non distributed parameters
546
547 Teuchos::ArrayRCP<const Scalar> p_data;
548 rcp_dynamic_cast<const Thyra::SpmdVectorBase<Scalar> >(paramVec,true)->getLocalData(Teuchos::ptrFromRef(p_data));
549
550 for (unsigned int j=0; j < parameters_[i]->scalar_value.size(); j++) {
551 parameters_[i]->scalar_value[j].baseValue = p_data[j];
552 parameters_[i]->scalar_value[j].family->setRealValueForAllTypes(parameters_[i]->scalar_value[j].baseValue);
553 }
554 }
555 else if ( paramVec!=Teuchos::null && parameters_[i]->is_distributed) {
556 // distributed parameters
557
558 std::string key = (*parameters_[i]->names)[0];
559 RCP<GlobalEvaluationData> ged = distrParamGlobalEvaluationData_.getDataObject(key);
560
561 TEUCHOS_ASSERT(ged!=Teuchos::null);
562
563 // cast to a LOCPair throwing an exception if the cast doesn't work.
564 RCP<LOCPair_GlobalEvaluationData> loc_pair_ged = rcp_dynamic_cast<LOCPair_GlobalEvaluationData>(ged);
565 RCP<ReadOnlyVector_GlobalEvaluationData> ro_ged = rcp_dynamic_cast<ReadOnlyVector_GlobalEvaluationData>(ged);
566 if(loc_pair_ged!=Teuchos::null) {
567 // cast to a ThyraObjContainer throwing an exception if the cast doesn't work.
568 RCP<ThyraObjContainer<Scalar> > th_ged = rcp_dynamic_cast<ThyraObjContainer<Scalar> >(loc_pair_ged->getGlobalLOC(),true);
569 th_ged->set_x_th(Teuchos::rcp_const_cast<Thyra::VectorBase<Scalar> >(paramVec));
570 }
571 else {
572 TEUCHOS_ASSERT(ro_ged!=Teuchos::null);
573 ro_ged->setOwnedVector(paramVec);
574 }
575 }
576 }
577
578 ae_inargs.addGlobalEvaluationData(distrParamGlobalEvaluationData_);
579
580 // here we are building a container, this operation is fast, simply allocating a struct
581 const RCP<panzer::ThyraObjContainer<Scalar> > thGlobalContainer =
582 Teuchos::rcp_dynamic_cast<panzer::ThyraObjContainer<Scalar> >(ae_inargs.container_);
583
584 TEUCHOS_ASSERT(!Teuchos::is_null(thGlobalContainer));
585
586 // Ghosted container objects are zeroed out below only if needed for
587 // a particular calculation. This makes it more efficient than
588 // zeroing out all objects in the container here.
589 // const RCP<panzer::ThyraObjContainer<Scalar> > thGhostedContainer =
590 // Teuchos::rcp_dynamic_cast<panzer::ThyraObjContainer<Scalar> >(ae_inargs.ghostedContainer_);
591
592 // Set the solution vector (currently all targets require solution).
593 // In the future we may move these into the individual cases below.
594 // A very subtle (and fragile) point: A non-null pointer in global
595 // container triggers export operations during fill. Also, the
596 // introduction of the container is forcing us to cast away const on
597 // arguments that should be const. Another reason to redesign
598 // LinearObjContainer layers.
599 thGlobalContainer->set_x_th(Teuchos::rcp_const_cast<Thyra::VectorBase<Scalar> >(x));
600 xContainer_->setOwnedVector(x);
601 ae_inargs.addGlobalEvaluationData("Solution Gather Container - X",xContainer_);
602
603 if (is_transient) {
604 thGlobalContainer->set_dxdt_th(Teuchos::rcp_const_cast<Thyra::VectorBase<Scalar> >(x_dot));
605 xdotContainer_->setOwnedVector(x_dot);
606 ae_inargs.addGlobalEvaluationData("Solution Gather Container - Xdot",xdotContainer_);
607 }
608
609 // Add tangent vectors for x and xdot to GlobalEvaluationData, one for each
610 // scalar parameter vector and parameter within that vector.
611 // Note: The keys for the global evaluation data containers for the tangent
612 // vectors are constructed in EquationSet_AddFieldDefaultImpl::
613 // buildAndRegisterGatherAndOrientationEvaluators().
614 int vIndex(0);
615 for (int i(0); i < num_param_vecs; ++i)
616 {
617 using std::string;
619 using Thyra::VectorBase;
621 if (not parameters_[i]->is_distributed)
622 {
623 auto dxdp = rcp_const_cast<VectorBase<Scalar>>
624 (inArgs.get_p(vIndex + num_param_vecs));
625 if (not dxdp.is_null())
626 {
627 // We need to cast away const because the object container requires
628 // non-const vectors.
629 auto dxdpBlock = rcp_dynamic_cast<ProductVectorBase<Scalar>>(dxdp);
630 int numParams(parameters_[i]->scalar_value.size());
631 for (int j(0); j < numParams; ++j)
632 {
633 RCP<ROVGED> dxdpContainer = lof_->buildReadOnlyDomainContainer();
634 dxdpContainer->setOwnedVector(dxdpBlock->getNonconstVectorBlock(j));
635 string name("X TANGENT GATHER CONTAINER: " +
636 (*parameters_[i]->names)[j]);
637 ae_inargs.addGlobalEvaluationData(name, dxdpContainer);
638 } // end loop over the parameters
639 } // end if (not dxdp.is_null())
640 if (build_transient_support_)
641 {
642 // We need to cast away const because the object container requires
643 // non-const vectors.
644 auto dxdotdp = rcp_const_cast<VectorBase<Scalar>>
645 (inArgs.get_p(vIndex + num_param_vecs + tangent_space_.size()));
646 if (not dxdotdp.is_null())
647 {
648 auto dxdotdpBlock =
649 rcp_dynamic_cast<ProductVectorBase<Scalar>>(dxdotdp);
650 int numParams(parameters_[i]->scalar_value.size());
651 for (int j(0); j < numParams; ++j)
652 {
653 RCP<ROVGED> dxdotdpContainer = lof_->buildReadOnlyDomainContainer();
654 dxdotdpContainer->setOwnedVector(
655 dxdotdpBlock->getNonconstVectorBlock(j));
656 string name("DXDT TANGENT GATHER CONTAINER: " +
657 (*parameters_[i]->names)[j]);
658 ae_inargs.addGlobalEvaluationData(name, dxdotdpContainer);
659 } // end loop over the parameters
660 } // end if (not dxdotdp.is_null())
661 } // end if (build_transient_support_)
662 ++vIndex;
663 } // end if (not parameters_[i]->is_distributed)
664 } // end loop over the parameter vectors
665} // end of setupAssemblyInArgs()
666
667// Private functions overridden from ModelEvaulatorDefaultBase
668
669
670template <typename Scalar>
671Thyra::ModelEvaluatorBase::OutArgs<Scalar>
673{
674 typedef Thyra::ModelEvaluatorBase MEB;
675
676 if(require_out_args_refresh_) {
677 MEB::OutArgsSetup<Scalar> outArgs;
678 outArgs.setModelEvalDescription(this->description());
679 outArgs.set_Np_Ng(num_me_parameters_, responses_.size());
680 outArgs.setSupports(MEB::OUT_ARG_f);
681 outArgs.setSupports(MEB::OUT_ARG_W_op);
682
683 // add in dg/dx (if appropriate)
684 for(std::size_t i=0;i<responses_.size();i++) {
685 typedef panzer::Traits::Jacobian RespEvalT;
686
687 // check dg/dx and add it in if appropriate
688 Teuchos::RCP<panzer::ResponseBase> respJacBase
689 = responseLibrary_->getResponse<RespEvalT>(responses_[i]->name);
690 if(respJacBase!=Teuchos::null) {
691 // cast is guranteed to succeed because of check in addResponse
692 Teuchos::RCP<panzer::ResponseMESupportBase<RespEvalT> > resp
693 = Teuchos::rcp_dynamic_cast<panzer::ResponseMESupportBase<RespEvalT> >(respJacBase);
694
695 // class must supppot a derivative
696 if(resp->supportsDerivative()) {
697 outArgs.setSupports(MEB::OUT_ARG_DgDx,i,MEB::DerivativeSupport(MEB::DERIV_MV_GRADIENT_FORM));
698
699
700 for(std::size_t p=0;p<parameters_.size();p++) {
701 if(parameters_[p]->is_distributed && parameters_[p]->global_indexer!=Teuchos::null)
702 outArgs.setSupports(MEB::OUT_ARG_DgDp,i,p,MEB::DerivativeSupport(MEB::DERIV_MV_GRADIENT_FORM));
703 if(!parameters_[p]->is_distributed)
704 outArgs.setSupports(MEB::OUT_ARG_DgDp,i,p,MEB::DerivativeSupport(MEB::DERIV_MV_JACOBIAN_FORM));
705 }
706 }
707 }
708 }
709
710 // setup parameter support
711 for(std::size_t p=0;p<parameters_.size();p++) {
712
713 if(!parameters_[p]->is_distributed)
714 outArgs.setSupports(MEB::OUT_ARG_DfDp,p,MEB::DerivativeSupport(MEB::DERIV_MV_BY_COL));
715 else if(parameters_[p]->is_distributed && parameters_[p]->global_indexer!=Teuchos::null)
716 outArgs.setSupports(MEB::OUT_ARG_DfDp,p,MEB::DerivativeSupport(MEB::DERIV_LINEAR_OP));
717 }
718
719 prototypeOutArgs_ = outArgs;
720 }
721
722 // we don't need to build it anymore
723 require_out_args_refresh_ = false;
724
725 return prototypeOutArgs_;
726}
727
728template <typename Scalar>
729Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
731create_W_op() const
732{
733 PANZER_FUNC_TIME_MONITOR("panzer::ModelEvaluator::create_W_op");
734 Teuchos::RCP<const ThyraObjFactory<Scalar> > tof
735 = Teuchos::rcp_dynamic_cast<const ThyraObjFactory<Scalar> >(lof_,true);
736
737 return tof->getThyraMatrix();
738}
739
740template <typename Scalar>
741Teuchos::RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> >
743get_W_factory() const
744{
745 return solverFactory_;
746}
747
748template <typename Scalar>
749Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
751create_DfDp_op(int p) const
752{
753 using Teuchos::RCP;
754 using Teuchos::rcp_dynamic_cast;
755
756 typedef Thyra::ModelEvaluatorBase MEB;
757
758 // The code below uses prototypeOutArgs_, so we need to make sure it is
759 // initialized before using it. This happens through createOutArgs(),
760 // however it may be allowable to call create_DfDp_op() before
761 // createOutArgs() is called. Thus we do this here if prototypeOutArgs_
762 // needs to be initialized.
763 //
764 // Previously this was handled in the TEUCHOS_ASSERT below through the call
765 // to Np(), however it isn't a good idea to include code in asserts that is
766 // required for proper execution (the asserts may be removed in an optimized
767 // build, for example).
768 if(require_out_args_refresh_) {
769 this->createOutArgs();
770 }
771
772 TEUCHOS_ASSERT(0<=p && p<Teuchos::as<int>(parameters_.size()));
773
774 // assert that DfDp is supported
775 const ParameterObject & po = *parameters_[p];
776
777 if(po.is_distributed && po.global_indexer!=Teuchos::null) {
778 TEUCHOS_ASSERT(prototypeOutArgs_.supports(MEB::OUT_ARG_DfDp,p).supports(MEB::DERIV_LINEAR_OP));
779
780 // for a distributed parameter, figure it out from the
781 // response library
782 RCP<Response_Residual<Traits::Jacobian> > response_jacobian
783 = rcp_dynamic_cast<Response_Residual<Traits::Jacobian> >(po.dfdp_rl->template getResponse<Traits::Jacobian>("RESIDUAL"));
784
785 return response_jacobian->allocateJacobian();
786 }
787 else if(!po.is_distributed) {
788 TEUCHOS_ASSERT(prototypeOutArgs_.supports(MEB::OUT_ARG_DfDp,p).supports(MEB::DERIV_MV_BY_COL));
789
790 // this is a scalar parameter (its easy to create!)
791 return Thyra::createMember(*get_f_space());
792 }
793
794 // shourld never get here
795 TEUCHOS_ASSERT(false);
796
797 return Teuchos::null;
798}
799
800template <typename Scalar>
802addParameter(const std::string & name,const Scalar & initialValue)
803{
804 Teuchos::Array<std::string> tmp_names;
805 tmp_names.push_back(name);
806
807 Teuchos::Array<Scalar> tmp_values;
808 tmp_values.push_back(initialValue);
809
810 return addParameter(tmp_names,tmp_values);
811}
812
813template <typename Scalar>
815addParameter(const Teuchos::Array<std::string> & names,
816 const Teuchos::Array<Scalar> & initialValues)
817{
818 using Teuchos::RCP;
819 using Teuchos::rcp;
820 using Teuchos::rcp_dynamic_cast;
821 using Teuchos::ptrFromRef;
822
823 TEUCHOS_ASSERT(names.size()==initialValues.size());
824
825 int parameter_index = parameters_.size();
826
827 // Create parameter object
828 RCP<ParameterObject> param = createScalarParameter(names,initialValues);
829 parameters_.push_back(param);
830
831 // Create vector space for parameter tangent vector
832 RCP< Thyra::VectorSpaceBase<double> > tan_space =
833 Thyra::multiVectorProductVectorSpace(x_space_, param->names->size());
834 tangent_space_.push_back(tan_space);
835
836 // The number of model evaluator parameters is the number of model parameters (parameters_.size()) plus a tangent
837 // vector for each scalar parameter (tangent_space_.size()) plus a tangent vector for xdot for each scalar parameter.
838 num_me_parameters_ += 2;
839 if (build_transient_support_)
840 ++num_me_parameters_;
841
842 require_in_args_refresh_ = true;
843 require_out_args_refresh_ = true;
844 this->resetDefaultBase();
845
846 return parameter_index;
847}
848
849template <typename Scalar>
851addDistributedParameter(const std::string & key,
852 const Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > & vs,
853 const Teuchos::RCP<GlobalEvaluationData> & ged,
854 const Teuchos::RCP<const Thyra::VectorBase<Scalar> > & initial,
855 const Teuchos::RCP<const GlobalIndexer> & ugi)
856{
857 distrParamGlobalEvaluationData_.addDataObject(key,ged);
858
859 int parameter_index = parameters_.size();
860 parameters_.push_back(createDistributedParameter(key,vs,initial,ugi));
861 ++num_me_parameters_;
862
863 require_in_args_refresh_ = true;
864 require_out_args_refresh_ = true;
865 this->resetDefaultBase();
866
867 return parameter_index;
868}
869
870template <typename Scalar>
872addNonParameterGlobalEvaluationData(const std::string & key,
873 const Teuchos::RCP<GlobalEvaluationData> & ged)
874{
875 nonParamGlobalEvaluationData_.addDataObject(key,ged);
876}
877
878template <typename Scalar>
880addFlexibleResponse(const std::string & responseName,
881 const std::vector<WorksetDescriptor> & wkst_desc,
882 const Teuchos::RCP<ResponseMESupportBuilderBase> & builder)
883{
884 // add a basic response, use x global indexer to define it
885 builder->setDerivativeInformation(lof_);
886
887 int respIndex = addResponse(responseName,wkst_desc,*builder);
888
889 // set the builder for this response
890 responses_[respIndex]->builder = builder;
891
892 return respIndex;
893}
894
895
896template <typename Scalar>
898applyDirichletBCs(const Teuchos::RCP<Thyra::VectorBase<Scalar> > & x,
899 const Teuchos::RCP<Thyra::VectorBase<Scalar> > & f) const
900{
901 using Teuchos::RCP;
902 using Teuchos::ArrayRCP;
903 using Teuchos::Array;
904 using Teuchos::tuple;
905 using Teuchos::rcp_dynamic_cast;
906
907 // if neccessary build a ghosted container
908 if(Teuchos::is_null(ghostedContainer_)) {
909 ghostedContainer_ = lof_->buildGhostedLinearObjContainer();
910 lof_->initializeGhostedContainer(panzer::LinearObjContainer::X |
911 panzer::LinearObjContainer::F,*ghostedContainer_);
912 }
913
915 ae_inargs.container_ = lof_->buildLinearObjContainer(); // we use a new global container
916 ae_inargs.ghostedContainer_ = ghostedContainer_; // we can reuse the ghosted container
917 ae_inargs.alpha = 0.0;
918 ae_inargs.beta = 1.0;
919 //TODO: is this really needed?
920 ae_inargs.step_size = 0.0;
921 ae_inargs.stage_number = 1.0;
922 ae_inargs.evaluate_transient_terms = false;
923 ae_inargs.addGlobalEvaluationData(nonParamGlobalEvaluationData_);
924 ae_inargs.addGlobalEvaluationData(distrParamGlobalEvaluationData_);
925
926 // this is the tempory target
927 lof_->initializeContainer(panzer::LinearObjContainer::F,*ae_inargs.container_);
928
929 // here we are building a container, this operation is fast, simply allocating a struct
930 const RCP<panzer::ThyraObjContainer<Scalar> > thGlobalContainer =
931 Teuchos::rcp_dynamic_cast<panzer::ThyraObjContainer<Scalar> >(ae_inargs.container_);
932
933 TEUCHOS_ASSERT(!Teuchos::is_null(thGlobalContainer));
934
935 // Ghosted container objects are zeroed out below only if needed for
936 // a particular calculation. This makes it more efficient than
937 // zeroing out all objects in the container here.
938 const RCP<panzer::ThyraObjContainer<Scalar> > thGhostedContainer =
939 Teuchos::rcp_dynamic_cast<panzer::ThyraObjContainer<Scalar> >(ae_inargs.ghostedContainer_);
940 Thyra::assign(thGhostedContainer->get_f_th().ptr(),0.0);
941
942 // Set the solution vector (currently all targets require solution).
943 // In the future we may move these into the individual cases below.
944 // A very subtle (and fragile) point: A non-null pointer in global
945 // container triggers export operations during fill. Also, the
946 // introduction of the container is forcing us to cast away const on
947 // arguments that should be const. Another reason to redesign
948 // LinearObjContainer layers.
949 thGlobalContainer->set_x_th(x);
950
951 // evaluate dirichlet boundary conditions
952 RCP<panzer::LinearObjContainer> counter
953 = ae_tm_.template getAsObject<panzer::Traits::Residual>()->evaluateOnlyDirichletBCs(ae_inargs);
954
955 // allocate the result container
956 RCP<panzer::LinearObjContainer> result = lof_->buildLinearObjContainer(); // we use a new global container
957
958 // stuff the evaluate boundary conditions into the f spot of the counter ... the x is already filled
959 Teuchos::rcp_dynamic_cast<panzer::ThyraObjContainer<Scalar> >(counter)->set_f_th(
960 thGlobalContainer->get_f_th());
961
962 // stuff the vector that needs applied dirichlet conditions in the the f spot of the result LOC
963 Teuchos::rcp_dynamic_cast<panzer::ThyraObjContainer<Scalar> >(result)->set_f_th(f);
964
965 // use the linear object factory to apply the result
966 lof_->applyDirichletBCs(*counter,*result);
967}
968
969template <typename Scalar>
971evalModel_D2gDx2(int respIndex,
972 const Thyra::ModelEvaluatorBase::InArgs<Scalar> & inArgs,
973 const Teuchos::RCP<const Thyra::VectorBase<Scalar> > & delta_x,
974 const Teuchos::RCP<Thyra::VectorBase<Scalar> > & D2gDx2) const
975{
976#ifdef Panzer_BUILD_HESSIAN_SUPPORT
977
978 // set model parameters from supplied inArgs
979 setParameters(inArgs);
980
981 {
982 std::string responseName = responses_[respIndex]->name;
983 Teuchos::RCP<panzer::ResponseMESupportBase<panzer::Traits::Hessian> > resp
984 = Teuchos::rcp_dynamic_cast<panzer::ResponseMESupportBase<panzer::Traits::Hessian> >(
985 responseLibrary_->getResponse<panzer::Traits::Hessian>(responseName));
986 resp->setDerivative(D2gDx2);
987 }
988
989 // setup all the assembly in arguments (this is parameters and
990 // x/x_dot). At this point with the exception of the one time dirichlet
991 // beta that is all thats neccessary.
993 setupAssemblyInArgs(inArgs,ae_inargs);
994
995 ae_inargs.beta = 1.0;
996
997 auto deltaXContainer = lof_->buildReadOnlyDomainContainer();
998 deltaXContainer->setOwnedVector(delta_x);
999 ae_inargs.addGlobalEvaluationData("DELTA_Solution Gather Container",deltaXContainer);
1000
1001 // evaluate responses
1002 responseLibrary_->addResponsesToInArgs<panzer::Traits::Hessian>(ae_inargs);
1003 responseLibrary_->evaluate<panzer::Traits::Hessian>(ae_inargs);
1004
1005 // reset parameters back to nominal values
1006 resetParameters();
1007#else
1008 (void)respIndex;
1009 (void)inArgs;
1010 (void)delta_x;
1011 (void)D2gDx2;
1012 TEUCHOS_ASSERT(false);
1013#endif
1014}
1015
1016template <typename Scalar>
1018evalModel_D2gDxDp(int respIndex,
1019 int pIndex,
1020 const Thyra::ModelEvaluatorBase::InArgs<Scalar> & inArgs,
1021 const Teuchos::RCP<const Thyra::VectorBase<Scalar> > & delta_p,
1022 const Teuchos::RCP<Thyra::VectorBase<Scalar> > & D2gDxDp) const
1023{
1024#ifdef Panzer_BUILD_HESSIAN_SUPPORT
1025
1026 // set model parameters from supplied inArgs
1027 setParameters(inArgs);
1028
1029 {
1030 std::string responseName = responses_[respIndex]->name;
1031 Teuchos::RCP<panzer::ResponseMESupportBase<panzer::Traits::Hessian> > resp
1032 = Teuchos::rcp_dynamic_cast<panzer::ResponseMESupportBase<panzer::Traits::Hessian> >(
1033 responseLibrary_->getResponse<panzer::Traits::Hessian>(responseName));
1034 resp->setDerivative(D2gDxDp);
1035 }
1036
1037 // setup all the assembly in arguments (this is parameters and
1038 // x/x_dot). At this point with the exception of the one time dirichlet
1039 // beta that is all thats neccessary.
1041 setupAssemblyInArgs(inArgs,ae_inargs);
1042
1043 ae_inargs.beta = 1.0;
1044 ae_inargs.second_sensitivities_name = (*parameters_[pIndex]->names)[0]; // distributed parameters have one name!
1045
1046 auto deltaPContainer = parameters_[pIndex]->dfdp_rl->getLinearObjFactory()->buildReadOnlyDomainContainer();
1047 deltaPContainer->setOwnedVector(delta_p);
1048 ae_inargs.addGlobalEvaluationData("DELTA_"+(*parameters_[pIndex]->names)[0],deltaPContainer);
1049
1050 // evaluate responses
1051 responseLibrary_->addResponsesToInArgs<panzer::Traits::Hessian>(ae_inargs);
1052 responseLibrary_->evaluate<panzer::Traits::Hessian>(ae_inargs);
1053
1054 // reset parameters back to nominal values
1055 resetParameters();
1056#else
1057 (void)respIndex;
1058 (void)pIndex;
1059 (void)inArgs;
1060 (void)delta_p;
1061 (void)D2gDxDp;
1062 TEUCHOS_ASSERT(false);
1063#endif
1064}
1065
1066template <typename Scalar>
1068evalModel_D2gDp2(int respIndex,
1069 int pIndex,
1070 const Thyra::ModelEvaluatorBase::InArgs<Scalar> & inArgs,
1071 const Teuchos::RCP<const Thyra::VectorBase<Scalar> > & delta_p,
1072 const Teuchos::RCP<Thyra::VectorBase<Scalar> > & D2gDp2) const
1073{
1074#ifdef Panzer_BUILD_HESSIAN_SUPPORT
1075
1076 // set model parameters from supplied inArgs
1077 setParameters(inArgs);
1078
1079 ResponseLibrary<Traits> & rLibrary = *parameters_[pIndex]->dgdp_rl;
1080
1081 {
1082 std::string responseName = responses_[respIndex]->name;
1083 Teuchos::RCP<panzer::ResponseMESupportBase<panzer::Traits::Hessian> > resp
1084 = Teuchos::rcp_dynamic_cast<panzer::ResponseMESupportBase<panzer::Traits::Hessian> >(
1085 rLibrary.getResponse<panzer::Traits::Hessian>(responseName));
1086 resp->setDerivative(D2gDp2);
1087 }
1088
1089 // setup all the assembly in arguments (this is parameters and
1090 // x/x_dot). At this point with the exception of the one time dirichlet
1091 // beta that is all thats neccessary.
1093 setupAssemblyInArgs(inArgs,ae_inargs);
1094
1095 ae_inargs.gather_seeds.push_back(1.0); // this assumes that gather point is always the zero index of
1096 // gather seeds
1097 ae_inargs.first_sensitivities_name = (*parameters_[pIndex]->names)[0]; // distributed parameters have one name!
1098 ae_inargs.second_sensitivities_name = (*parameters_[pIndex]->names)[0]; // distributed parameters have one name!
1099
1100 auto deltaPContainer = parameters_[pIndex]->dfdp_rl->getLinearObjFactory()->buildReadOnlyDomainContainer();
1101 deltaPContainer->setOwnedVector(delta_p);
1102 ae_inargs.addGlobalEvaluationData("DELTA_"+(*parameters_[pIndex]->names)[0],deltaPContainer);
1103
1104 // evaluate responses
1105 rLibrary.addResponsesToInArgs<panzer::Traits::Hessian>(ae_inargs);
1106 rLibrary.evaluate<panzer::Traits::Hessian>(ae_inargs);
1107
1108 // reset parameters back to nominal values
1109 resetParameters();
1110#else
1111 (void)respIndex;
1112 (void)pIndex;
1113 (void)inArgs;
1114 (void)delta_p;
1115 (void)D2gDp2;
1116 TEUCHOS_ASSERT(false);
1117#endif
1118}
1119
1120template <typename Scalar>
1122evalModel_D2gDpDx(int respIndex,
1123 int pIndex,
1124 const Thyra::ModelEvaluatorBase::InArgs<Scalar> & inArgs,
1125 const Teuchos::RCP<const Thyra::VectorBase<Scalar> > & delta_x,
1126 const Teuchos::RCP<Thyra::VectorBase<Scalar> > & D2gDpDx) const
1127{
1128#ifdef Panzer_BUILD_HESSIAN_SUPPORT
1129
1130 // set model parameters from supplied inArgs
1131 setParameters(inArgs);
1132
1133 ResponseLibrary<Traits> & rLibrary = *parameters_[pIndex]->dgdp_rl;
1134
1135 {
1136 std::string responseName = responses_[respIndex]->name;
1137 Teuchos::RCP<panzer::ResponseMESupportBase<panzer::Traits::Hessian> > resp
1138 = Teuchos::rcp_dynamic_cast<panzer::ResponseMESupportBase<panzer::Traits::Hessian> >(
1139 rLibrary.getResponse<panzer::Traits::Hessian>(responseName));
1140 resp->setDerivative(D2gDpDx);
1141 }
1142
1143 // setup all the assembly in arguments (this is parameters and
1144 // x/x_dot). At this point with the exception of the one time dirichlet
1145 // beta that is all thats neccessary.
1147 setupAssemblyInArgs(inArgs,ae_inargs);
1148
1149 ae_inargs.gather_seeds.push_back(1.0); // this assumes that gather point is always the zero index of
1150 // gather seeds
1151 ae_inargs.first_sensitivities_name = (*parameters_[pIndex]->names)[0]; // distributed parameters have one name!
1152 ae_inargs.second_sensitivities_name = "";
1153
1154 auto deltaXContainer = lof_->buildReadOnlyDomainContainer();
1155 deltaXContainer->setOwnedVector(delta_x);
1156 ae_inargs.addGlobalEvaluationData("DELTA_Solution Gather Container",deltaXContainer);
1157
1158 // evaluate responses
1159 rLibrary.addResponsesToInArgs<panzer::Traits::Hessian>(ae_inargs);
1160 rLibrary.evaluate<panzer::Traits::Hessian>(ae_inargs);
1161
1162 // reset parameters back to nominal values
1163 resetParameters();
1164#else
1165 (void)respIndex;
1166 (void)pIndex;
1167 (void)inArgs;
1168 (void)delta_x;
1169 (void)D2gDpDx;
1170 TEUCHOS_ASSERT(false);
1171#endif
1172}
1173
1174template <typename Scalar>
1176evalModel_D2fDx2(const Thyra::ModelEvaluatorBase::InArgs<Scalar> & inArgs,
1177 const Teuchos::RCP<const Thyra::VectorBase<Scalar> > & delta_x,
1178 const Teuchos::RCP<Thyra::LinearOpBase<Scalar> > & D2fDx2) const
1179{
1180#ifdef Panzer_BUILD_HESSIAN_SUPPORT
1181
1182 using Teuchos::RCP;
1183 using Teuchos::ArrayRCP;
1184 using Teuchos::Array;
1185 using Teuchos::tuple;
1186 using Teuchos::rcp_dynamic_cast;
1187
1188 typedef Thyra::ModelEvaluatorBase MEB;
1189
1190 // Transient or steady-state evaluation is determined by the x_dot
1191 // vector. If this RCP is null, then we are doing a steady-state
1192 // fill.
1193 bool is_transient = false;
1194 if (inArgs.supports(MEB::IN_ARG_x_dot ))
1195 is_transient = !Teuchos::is_null(inArgs.get_x_dot());
1196
1197 // Make sure construction built in transient support
1198 TEUCHOS_TEST_FOR_EXCEPTION(is_transient && !build_transient_support_, std::runtime_error,
1199 "ModelEvaluator was not built with transient support enabled!");
1200
1201 //
1202 // Get the output arguments
1203 //
1204 const RCP<Thyra::LinearOpBase<Scalar> > W_out = D2fDx2;
1205
1206 // setup all the assembly in arguments (this is parameters and
1207 // x/x_dot). At this point with the exception of the one time dirichlet
1208 // beta that is all thats neccessary.
1210 setupAssemblyInArgs(inArgs,ae_inargs);
1211
1212 auto deltaXContainer = lof_->buildReadOnlyDomainContainer();
1213 deltaXContainer->setOwnedVector(delta_x);
1214 ae_inargs.addGlobalEvaluationData("DELTA_Solution Gather Container",deltaXContainer);
1215
1216 // set model parameters from supplied inArgs
1217 setParameters(inArgs);
1218
1219 // handle application of the one time dirichlet beta in the
1220 // assembly engine. Note that this has to be set explicitly
1221 // each time because this badly breaks encapsulation. Essentially
1222 // we must work around the model evaluator abstraction!
1223 if(oneTimeDirichletBeta_on_) {
1224 ae_inargs.dirichlet_beta = oneTimeDirichletBeta_;
1225 ae_inargs.apply_dirichlet_beta = true;
1226
1227 oneTimeDirichletBeta_on_ = false;
1228 }
1229
1230 // here we are building a container, this operation is fast, simply allocating a struct
1231 const RCP<panzer::ThyraObjContainer<Scalar> > thGlobalContainer =
1232 Teuchos::rcp_dynamic_cast<panzer::ThyraObjContainer<Scalar> >(ae_inargs.container_);
1233 const RCP<panzer::ThyraObjContainer<Scalar> > thGhostedContainer =
1234 Teuchos::rcp_dynamic_cast<panzer::ThyraObjContainer<Scalar> >(ae_inargs.ghostedContainer_);
1235
1236 {
1237 PANZER_FUNC_TIME_MONITOR("panzer::ModelEvaluator::evalModel(D2fDx2)");
1238
1239 // this dummy nonsense is needed only for scattering dirichlet conditions
1240 RCP<Thyra::VectorBase<Scalar> > dummy_f = Thyra::createMember(f_space_);
1241 thGlobalContainer->set_f_th(dummy_f);
1242 thGlobalContainer->set_A_th(W_out);
1243
1244 // Zero values in ghosted container objects
1245 thGhostedContainer->initializeMatrix(0.0);
1246
1247 ae_tm_.template getAsObject<panzer::Traits::Hessian>()->evaluate(ae_inargs);
1248 }
1249
1250 // HACK: set A to null before calling responses to avoid touching the
1251 // the Jacobian after it has been properly assembled. Should be fixed
1252 // by using a modified version of ae_inargs instead.
1253 thGlobalContainer->set_A_th(Teuchos::null);
1254
1255 // Holding a rcp to f produces a seg fault in Rythmos when the next
1256 // f comes in and the resulting dtor is called. Need to discuss
1257 // with Ross. Clearing all references here works!
1258
1259 thGlobalContainer->set_x_th(Teuchos::null);
1260 thGlobalContainer->set_dxdt_th(Teuchos::null);
1261 thGlobalContainer->set_f_th(Teuchos::null);
1262 thGlobalContainer->set_A_th(Teuchos::null);
1263
1264 // reset parameters back to nominal values
1265 resetParameters();
1266#else
1267 (void)inArgs;
1268 (void)delta_x;
1269 (void)D2fDx2;
1270 TEUCHOS_ASSERT(false);
1271#endif
1272}
1273
1274template <typename Scalar>
1276evalModel_D2fDxDp(int pIndex,
1277 const Thyra::ModelEvaluatorBase::InArgs<Scalar> & inArgs,
1278 const Teuchos::RCP<const Thyra::VectorBase<Scalar> > & delta_p,
1279 const Teuchos::RCP<Thyra::LinearOpBase<Scalar> > & D2fDxDp) const
1280{
1281#ifdef Panzer_BUILD_HESSIAN_SUPPORT
1282
1283 using Teuchos::RCP;
1284 using Teuchos::ArrayRCP;
1285 using Teuchos::Array;
1286 using Teuchos::tuple;
1287 using Teuchos::rcp_dynamic_cast;
1288
1289 typedef Thyra::ModelEvaluatorBase MEB;
1290
1291 // Transient or steady-state evaluation is determined by the x_dot
1292 // vector. If this RCP is null, then we are doing a steady-state
1293 // fill.
1294 bool is_transient = false;
1295 if (inArgs.supports(MEB::IN_ARG_x_dot ))
1296 is_transient = !Teuchos::is_null(inArgs.get_x_dot());
1297
1298 // Make sure construction built in transient support
1299 TEUCHOS_TEST_FOR_EXCEPTION(is_transient && !build_transient_support_, std::runtime_error,
1300 "ModelEvaluator was not built with transient support enabled!");
1301
1302 //
1303 // Get the output arguments
1304 //
1305 const RCP<Thyra::LinearOpBase<Scalar> > W_out = D2fDxDp;
1306
1307 // setup all the assembly in arguments (this is parameters and
1308 // x/x_dot). At this point with the exception of the one time dirichlet
1309 // beta that is all thats neccessary.
1311 setupAssemblyInArgs(inArgs,ae_inargs);
1312
1313 ae_inargs.second_sensitivities_name = (*parameters_[pIndex]->names)[0]; // distributed parameters have one name!
1314
1315 auto deltaPContainer = parameters_[pIndex]->dfdp_rl->getLinearObjFactory()->buildReadOnlyDomainContainer();
1316 deltaPContainer->setOwnedVector(delta_p);
1317 ae_inargs.addGlobalEvaluationData("DELTA_"+(*parameters_[pIndex]->names)[0],deltaPContainer);
1318
1319 // set model parameters from supplied inArgs
1320 setParameters(inArgs);
1321
1322 // handle application of the one time dirichlet beta in the
1323 // assembly engine. Note that this has to be set explicitly
1324 // each time because this badly breaks encapsulation. Essentially
1325 // we must work around the model evaluator abstraction!
1326 if(oneTimeDirichletBeta_on_) {
1327 ae_inargs.dirichlet_beta = oneTimeDirichletBeta_;
1328 ae_inargs.apply_dirichlet_beta = true;
1329
1330 oneTimeDirichletBeta_on_ = false;
1331 }
1332
1333 // here we are building a container, this operation is fast, simply allocating a struct
1334 const RCP<panzer::ThyraObjContainer<Scalar> > thGlobalContainer =
1335 Teuchos::rcp_dynamic_cast<panzer::ThyraObjContainer<Scalar> >(ae_inargs.container_);
1336 const RCP<panzer::ThyraObjContainer<Scalar> > thGhostedContainer =
1337 Teuchos::rcp_dynamic_cast<panzer::ThyraObjContainer<Scalar> >(ae_inargs.ghostedContainer_);
1338
1339 {
1340 PANZER_FUNC_TIME_MONITOR("panzer::ModelEvaluator::evalModel(D2fDxDp)");
1341
1342 // this dummy nonsense is needed only for scattering dirichlet conditions
1343 RCP<Thyra::VectorBase<Scalar> > dummy_f = Thyra::createMember(f_space_);
1344 thGlobalContainer->set_f_th(dummy_f);
1345 thGlobalContainer->set_A_th(W_out);
1346
1347 // Zero values in ghosted container objects
1348 thGhostedContainer->initializeMatrix(0.0);
1349
1350 ae_tm_.template getAsObject<panzer::Traits::Hessian>()->evaluate(ae_inargs);
1351 }
1352
1353 // HACK: set A to null before calling responses to avoid touching the
1354 // the Jacobian after it has been properly assembled. Should be fixed
1355 // by using a modified version of ae_inargs instead.
1356 thGlobalContainer->set_A_th(Teuchos::null);
1357
1358 // Holding a rcp to f produces a seg fault in Rythmos when the next
1359 // f comes in and the resulting dtor is called. Need to discuss
1360 // with Ross. Clearing all references here works!
1361
1362 thGlobalContainer->set_x_th(Teuchos::null);
1363 thGlobalContainer->set_dxdt_th(Teuchos::null);
1364 thGlobalContainer->set_f_th(Teuchos::null);
1365 thGlobalContainer->set_A_th(Teuchos::null);
1366
1367 // reset parameters back to nominal values
1368 resetParameters();
1369#else
1370 (void)pIndex;
1371 (void)inArgs;
1372 (void)delta_p;
1373 (void)D2fDxDp;
1374 TEUCHOS_ASSERT(false);
1375#endif
1376}
1377
1378template <typename Scalar>
1380evalModel_D2fDpDx(int pIndex,
1381 const Thyra::ModelEvaluatorBase::InArgs<Scalar> & inArgs,
1382 const Teuchos::RCP<const Thyra::VectorBase<Scalar> > & delta_x,
1383 const Teuchos::RCP<Thyra::LinearOpBase<Scalar> > & D2fDpDx) const
1384{
1385#ifdef Panzer_BUILD_HESSIAN_SUPPORT
1386 using Teuchos::RCP;
1387 using Teuchos::rcp_dynamic_cast;
1388 using Teuchos::null;
1389
1390 // parameter is not distributed
1391 TEUCHOS_ASSERT(parameters_[pIndex]->is_distributed);
1392
1393 // parameter is distributed but has no global indexer.
1394 // thus the user doesn't want sensitivities!
1395 TEUCHOS_ASSERT(parameters_[pIndex]->dfdp_rl!=null);
1396
1397 ResponseLibrary<Traits> & rLibrary = *parameters_[pIndex]->dfdp_rl;
1398
1399 // get the response and tell it to fill the derivative operator
1400 RCP<Response_Residual<Traits::Hessian> > response_hessian =
1401 rcp_dynamic_cast<Response_Residual<Traits::Hessian> >(rLibrary.getResponse<Traits::Hessian>("RESIDUAL"));
1402 response_hessian->setHessian(D2fDpDx);
1403
1404 // setup all the assembly in arguments (this is parameters and x/x_dot).
1405 // make sure the correct seeding is performed
1407 setupAssemblyInArgs(inArgs,ae_inargs);
1408
1409 auto deltaXContainer = lof_->buildReadOnlyDomainContainer();
1410 deltaXContainer->setOwnedVector(delta_x);
1411 ae_inargs.addGlobalEvaluationData("DELTA_Solution Gather Container",deltaXContainer);
1412
1413 ae_inargs.gather_seeds.push_back(1.0); // this assumes that gather point is always the zero index of
1414 // gather seeds
1415 ae_inargs.first_sensitivities_name = (*parameters_[pIndex]->names)[0]; // distributed parameters have one name!
1416 ae_inargs.second_sensitivities_name = "";
1417
1418 rLibrary.addResponsesToInArgs<Traits::Hessian>(ae_inargs);
1419 rLibrary.evaluate<Traits::Hessian>(ae_inargs);
1420#else
1421 (void)pIndex;
1422 (void)inArgs;
1423 (void)delta_x;
1424 (void)D2fDpDx;
1425 TEUCHOS_ASSERT(false);
1426#endif
1427}
1428
1429template <typename Scalar>
1431evalModel_D2fDp2(int pIndex,
1432 const Thyra::ModelEvaluatorBase::InArgs<Scalar> & inArgs,
1433 const Teuchos::RCP<const Thyra::VectorBase<Scalar> > & delta_p,
1434 const Teuchos::RCP<Thyra::LinearOpBase<Scalar> > & D2fDp2) const
1435{
1436#ifdef Panzer_BUILD_HESSIAN_SUPPORT
1437 using Teuchos::RCP;
1438 using Teuchos::rcp_dynamic_cast;
1439 using Teuchos::null;
1440
1441 // parameter is not distributed
1442 TEUCHOS_ASSERT(parameters_[pIndex]->is_distributed);
1443
1444 // parameter is distributed but has no global indexer.
1445 // thus the user doesn't want sensitivities!
1446 TEUCHOS_ASSERT(parameters_[pIndex]->dfdp_rl!=null);
1447
1448 ResponseLibrary<Traits> & rLibrary = *parameters_[pIndex]->dfdp_rl;
1449
1450 // get the response and tell it to fill the derivative operator
1451 RCP<Response_Residual<Traits::Hessian> > response_hessian =
1452 rcp_dynamic_cast<Response_Residual<Traits::Hessian> >(rLibrary.getResponse<Traits::Hessian>("RESIDUAL"));
1453 response_hessian->setHessian(D2fDp2);
1454
1455 // setup all the assembly in arguments (this is parameters and x/x_dot).
1456 // make sure the correct seeding is performed
1458 setupAssemblyInArgs(inArgs,ae_inargs);
1459
1460 auto deltaPContainer = parameters_[pIndex]->dfdp_rl->getLinearObjFactory()->buildReadOnlyDomainContainer();
1461 deltaPContainer->setOwnedVector(delta_p);
1462 ae_inargs.addGlobalEvaluationData("DELTA_"+(*parameters_[pIndex]->names)[0],deltaPContainer);
1463
1464 ae_inargs.gather_seeds.push_back(1.0); // this assumes that gather point is always the zero index of
1465 // gather seeds
1466 ae_inargs.first_sensitivities_name = (*parameters_[pIndex]->names)[0]; // distributed parameters have one name!
1467 ae_inargs.second_sensitivities_name = (*parameters_[pIndex]->names)[0]; // distributed parameters have one name!
1468
1469 rLibrary.addResponsesToInArgs<Traits::Hessian>(ae_inargs);
1470 rLibrary.evaluate<Traits::Hessian>(ae_inargs);
1471#else
1472 (void)pIndex;
1473 (void)inArgs;
1474 (void)delta_p;
1475 (void)D2fDp2;
1476 TEUCHOS_ASSERT(false);
1477#endif
1478}
1479
1480template <typename Scalar>
1482evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
1483 const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
1484{
1485 evalModelImpl_basic(inArgs,outArgs);
1486
1487 // evaluate responses...uses the stored assembly arguments and containers
1488 if(required_basic_g(outArgs))
1489 evalModelImpl_basic_g(inArgs,outArgs);
1490
1491 // evaluate response derivatives
1492 if(required_basic_dgdx(outArgs))
1493 evalModelImpl_basic_dgdx(inArgs,outArgs);
1494
1495 // evaluate response derivatives to scalar parameters
1496 if(required_basic_dgdp_scalar(outArgs))
1497 evalModelImpl_basic_dgdp_scalar(inArgs,outArgs);
1498
1499 // evaluate response derivatives to distributed parameters
1500 if(required_basic_dgdp_distro(outArgs))
1501 evalModelImpl_basic_dgdp_distro(inArgs,outArgs);
1502
1503 if(required_basic_dfdp_scalar(outArgs)) {
1504 if (do_fd_dfdp_)
1505 evalModelImpl_basic_dfdp_scalar_fd(inArgs,outArgs);
1506 else
1507 evalModelImpl_basic_dfdp_scalar(inArgs,outArgs);
1508 }
1509
1510 if(required_basic_dfdp_distro(outArgs))
1511 evalModelImpl_basic_dfdp_distro(inArgs,outArgs);
1512}
1513
1514template <typename Scalar>
1516evalModelImpl_basic(const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
1517 const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
1518{
1519 using Teuchos::RCP;
1520 using Teuchos::ArrayRCP;
1521 using Teuchos::Array;
1522 using Teuchos::tuple;
1523 using Teuchos::rcp_dynamic_cast;
1524
1525 typedef Thyra::ModelEvaluatorBase MEB;
1526
1527 // Transient or steady-state evaluation is determined by the x_dot
1528 // vector. If this RCP is null, then we are doing a steady-state
1529 // fill.
1530 bool is_transient = false;
1531 if (inArgs.supports(MEB::IN_ARG_x_dot ))
1532 is_transient = !Teuchos::is_null(inArgs.get_x_dot());
1533
1534 // Make sure construction built in transient support
1535 TEUCHOS_TEST_FOR_EXCEPTION(is_transient && !build_transient_support_, std::runtime_error,
1536 "ModelEvaluator was not built with transient support enabled!");
1537
1538 //
1539 // Get the output arguments
1540 //
1541 const RCP<Thyra::VectorBase<Scalar> > f_out = outArgs.get_f();
1542 const RCP<Thyra::LinearOpBase<Scalar> > W_out = outArgs.get_W_op();
1543
1544 // see if the user wants us to do anything
1545 if(Teuchos::is_null(f_out) && Teuchos::is_null(W_out) ) {
1546 return;
1547 }
1548
1549 // setup all the assembly in arguments (this is parameters and
1550 // x/x_dot). At this point with the exception of the one time dirichlet
1551 // beta that is all thats neccessary.
1553 setupAssemblyInArgs(inArgs,ae_inargs);
1554
1555 // set model parameters from supplied inArgs
1556 setParameters(inArgs);
1557
1558 // handle application of the one time dirichlet beta in the
1559 // assembly engine. Note that this has to be set explicitly
1560 // each time because this badly breaks encapsulation. Essentially
1561 // we must work around the model evaluator abstraction!
1562 if(oneTimeDirichletBeta_on_) {
1563 ae_inargs.dirichlet_beta = oneTimeDirichletBeta_;
1564 ae_inargs.apply_dirichlet_beta = true;
1565
1566 oneTimeDirichletBeta_on_ = false;
1567 }
1568
1569 // here we are building a container, this operation is fast, simply allocating a struct
1570 const RCP<panzer::ThyraObjContainer<Scalar> > thGlobalContainer =
1571 Teuchos::rcp_dynamic_cast<panzer::ThyraObjContainer<Scalar> >(ae_inargs.container_);
1572 const RCP<panzer::ThyraObjContainer<Scalar> > thGhostedContainer =
1573 Teuchos::rcp_dynamic_cast<panzer::ThyraObjContainer<Scalar> >(ae_inargs.ghostedContainer_);
1574
1575 if (!Teuchos::is_null(f_out) && !Teuchos::is_null(W_out)) {
1576 PANZER_FUNC_TIME_MONITOR("panzer::ModelEvaluator::evalModel(f and J)");
1577
1578 // only add auxiliary global data if Jacobian is being formed
1579 ae_inargs.addGlobalEvaluationData(nonParamGlobalEvaluationData_);
1580
1581 // Set the targets
1582 thGlobalContainer->set_f_th(f_out);
1583 thGlobalContainer->set_A_th(W_out);
1584
1585 // Zero values in ghosted container objects
1586 Thyra::assign(thGhostedContainer->get_f_th().ptr(),0.0);
1587 thGhostedContainer->initializeMatrix(0.0);
1588
1589 ae_tm_.template getAsObject<panzer::Traits::Jacobian>()->evaluate(ae_inargs);
1590 }
1591 else if(!Teuchos::is_null(f_out) && Teuchos::is_null(W_out)) {
1592
1593 PANZER_FUNC_TIME_MONITOR("panzer::ModelEvaluator::evalModel(f)");
1594
1595 // don't add auxiliary global data if Jacobian is not computed.
1596 // this leads to zeroing of aux ops in special cases.
1597
1598 thGlobalContainer->set_f_th(f_out);
1599
1600 // Zero values in ghosted container objects
1601 Thyra::assign(thGhostedContainer->get_f_th().ptr(),0.0);
1602
1603 ae_tm_.template getAsObject<panzer::Traits::Residual>()->evaluate(ae_inargs);
1604 }
1605 else if(Teuchos::is_null(f_out) && !Teuchos::is_null(W_out)) {
1606
1607 PANZER_FUNC_TIME_MONITOR("panzer::ModelEvaluator::evalModel(J)");
1608
1609 // only add auxiliary global data if Jacobian is being formed
1610 ae_inargs.addGlobalEvaluationData(nonParamGlobalEvaluationData_);
1611
1612 // this dummy nonsense is needed only for scattering dirichlet conditions
1613 RCP<Thyra::VectorBase<Scalar> > dummy_f = Thyra::createMember(f_space_);
1614 thGlobalContainer->set_f_th(dummy_f);
1615 thGlobalContainer->set_A_th(W_out);
1616
1617 // Zero values in ghosted container objects
1618 thGhostedContainer->initializeMatrix(0.0);
1619
1620 ae_tm_.template getAsObject<panzer::Traits::Jacobian>()->evaluate(ae_inargs);
1621 }
1622
1623 // HACK: set A to null before calling responses to avoid touching the
1624 // the Jacobian after it has been properly assembled. Should be fixed
1625 // by using a modified version of ae_inargs instead.
1626 thGlobalContainer->set_A_th(Teuchos::null);
1627
1628 // Holding a rcp to f produces a seg fault in Rythmos when the next
1629 // f comes in and the resulting dtor is called. Need to discuss
1630 // with Ross. Clearing all references here works!
1631
1632 thGlobalContainer->set_x_th(Teuchos::null);
1633 thGlobalContainer->set_dxdt_th(Teuchos::null);
1634 thGlobalContainer->set_f_th(Teuchos::null);
1635 thGlobalContainer->set_A_th(Teuchos::null);
1636
1637 // reset parameters back to nominal values
1638 resetParameters();
1639
1640 const bool writeToFile = false;
1641 if (writeToFile && nonnull(W_out)) {
1642 const auto check_blocked = Teuchos::rcp_dynamic_cast<::Thyra::BlockedLinearOpBase<double> >(W_out,false);
1643 if (check_blocked) {
1644 const int numBlocks = check_blocked->productDomain()->numBlocks();
1645 const int rangeBlocks = check_blocked->productRange()->numBlocks();
1646 TEUCHOS_ASSERT(numBlocks == rangeBlocks); // not true for optimization
1647 for (int row=0; row < numBlocks; ++row) {
1648 for (int col=0; col < numBlocks; ++col) {
1649 using LO = panzer::LocalOrdinal;
1650 using GO = panzer::GlobalOrdinal;
1651 using NodeT = panzer::TpetraNodeType;
1652 const auto thyraTpetraOperator = Teuchos::rcp_dynamic_cast<::Thyra::TpetraLinearOp<double,LO,GO,NodeT>>(check_blocked->getNonconstBlock(row,col),true);
1653 const auto tpetraCrsMatrix = Teuchos::rcp_dynamic_cast<Tpetra::CrsMatrix<double,LO,GO,NodeT>>(thyraTpetraOperator->getTpetraOperator(),true);
1654 tpetraCrsMatrix->print(std::cout);
1655 std::stringstream ss;
1656 ss << "W_out_" << write_matrix_count_ << ".rank_" << tpetraCrsMatrix->getMap()->getComm()->getRank() << ".block_" << row << "_" << col << ".txt";
1657 std::fstream fs(ss.str().c_str(),std::fstream::out|std::fstream::trunc);
1658 Teuchos::FancyOStream fos(Teuchos::rcpFromRef(fs));
1659 tpetraCrsMatrix->describe(fos,Teuchos::VERB_EXTREME);
1660 fs.close();
1661 }
1662 }
1663 }
1664 else {
1665 using LO = panzer::LocalOrdinal;
1666 using GO = panzer::GlobalOrdinal;
1667 using NodeT = panzer::TpetraNodeType;
1668 const auto thyraTpetraOperator = Teuchos::rcp_dynamic_cast<::Thyra::TpetraLinearOp<double,LO,GO,NodeT>>(W_out,true);
1669 const auto tpetraCrsMatrix = Teuchos::rcp_dynamic_cast<Tpetra::CrsMatrix<double,LO,GO,NodeT>>(thyraTpetraOperator->getTpetraOperator(),true);
1670 tpetraCrsMatrix->print(std::cout);
1671 std::stringstream ss;
1672 ss << "W_out_" << write_matrix_count_ << ".rank_" << tpetraCrsMatrix->getMap()->getComm()->getRank() << ".txt";
1673 std::fstream fs(ss.str().c_str(),std::fstream::out|std::fstream::trunc);
1674 Teuchos::FancyOStream fos(Teuchos::rcpFromRef(fs));
1675 tpetraCrsMatrix->describe(fos,Teuchos::VERB_EXTREME);
1676 fs.close();
1677 }
1678 ++write_matrix_count_;
1679 }
1680
1681}
1682
1683template <typename Scalar>
1685evalModelImpl_basic_g(const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
1686 const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
1687{
1688 // optional sanity check
1689 // TEUCHOS_ASSERT(required_basic_g(outArgs));
1690
1691 // setup all the assembly in arguments (this is parameters and
1692 // x/x_dot). At this point with the exception of the one time dirichlet
1693 // beta that is all thats neccessary.
1695 setupAssemblyInArgs(inArgs,ae_inargs);
1696
1697 // set model parameters from supplied inArgs
1698 setParameters(inArgs);
1699
1700 for(std::size_t i=0;i<responses_.size();i++) {
1701 Teuchos::RCP<Thyra::VectorBase<Scalar> > vec = outArgs.get_g(i);
1702 if(vec!=Teuchos::null) {
1703 std::string responseName = responses_[i]->name;
1704 Teuchos::RCP<panzer::ResponseMESupportBase<panzer::Traits::Residual> > resp
1705 = Teuchos::rcp_dynamic_cast<panzer::ResponseMESupportBase<panzer::Traits::Residual> >(
1706 responseLibrary_->getResponse<panzer::Traits::Residual>(responseName));
1707 resp->setVector(vec);
1708 }
1709 }
1710
1711 // evaluator responses
1712 responseLibrary_->addResponsesToInArgs<panzer::Traits::Residual>(ae_inargs);
1713 responseLibrary_->evaluate<panzer::Traits::Residual>(ae_inargs);
1714
1715 // reset parameters back to nominal values
1716 resetParameters();
1717}
1718
1719template <typename Scalar>
1720void
1722evalModelImpl_basic_dgdx(const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
1723 const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
1724{
1725 typedef Thyra::ModelEvaluatorBase MEB;
1726
1727 // optional sanity check
1728 TEUCHOS_ASSERT(required_basic_dgdx(outArgs));
1729
1730 // set model parameters from supplied inArgs
1731 setParameters(inArgs);
1732
1733 for(std::size_t i=0;i<responses_.size();i++) {
1734 // get "Vector" out of derivative, if its something else, throw an exception
1735 MEB::Derivative<Scalar> deriv = outArgs.get_DgDx(i);
1736 if(deriv.isEmpty())
1737 continue;
1738
1739 Teuchos::RCP<Thyra::MultiVectorBase<Scalar> > vec = deriv.getMultiVector();
1740
1741 if(vec!=Teuchos::null) {
1742
1743 std::string responseName = responses_[i]->name;
1744 Teuchos::RCP<panzer::ResponseMESupportBase<panzer::Traits::Jacobian> > resp
1745 = Teuchos::rcp_dynamic_cast<panzer::ResponseMESupportBase<panzer::Traits::Jacobian> >(
1746 responseLibrary_->getResponse<panzer::Traits::Jacobian>(responseName));
1747 resp->setDerivative(vec);
1748 }
1749 }
1750
1751 // setup all the assembly in arguments (this is parameters and
1752 // x/x_dot). At this point with the exception of the one time dirichlet
1753 // beta that is all thats neccessary.
1755 setupAssemblyInArgs(inArgs,ae_inargs);
1756
1757 // evaluate responses
1758 responseLibrary_->addResponsesToInArgs<panzer::Traits::Jacobian>(ae_inargs);
1759 responseLibrary_->evaluate<panzer::Traits::Jacobian>(ae_inargs);
1760
1761 // reset parameters back to nominal values
1762 resetParameters();
1763}
1764
1765template <typename Scalar>
1766void
1768evalModelImpl_basic_dgdp_scalar(const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
1769 const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
1770{
1771 using Teuchos::RCP;
1772 using Teuchos::rcp;
1773 using Teuchos::rcp_dynamic_cast;
1774
1775 typedef Thyra::ModelEvaluatorBase MEB;
1776
1777 // optional sanity check
1778 TEUCHOS_ASSERT(required_basic_dgdp_scalar(outArgs));
1779
1780 // First find all of the active parameters for all responses
1781 std::vector<std::string> activeParameterNames;
1782 std::vector<int> activeParameters;
1783 int totalParameterCount = 0;
1784 for(std::size_t j=0; j<parameters_.size(); j++) {
1785
1786 // skip non-scalar parameters
1787 if(parameters_[j]->is_distributed)
1788 continue;
1789
1790 bool is_active = false;
1791 for(std::size_t i=0;i<responses_.size(); i++) {
1792
1793 MEB::Derivative<Scalar> deriv = outArgs.get_DgDp(i,j);
1794 if(deriv.isEmpty())
1795 continue;
1796
1797 Teuchos::RCP<Thyra::MultiVectorBase<Scalar> > vec = deriv.getMultiVector();
1798 if(vec!=Teuchos::null) {
1799 // get the response and tell it to fill the derivative vector
1800 std::string responseName = responses_[i]->name;
1801 RCP<panzer::ResponseMESupportBase<panzer::Traits::Tangent> > resp =
1802 rcp_dynamic_cast<panzer::ResponseMESupportBase<panzer::Traits::Tangent> >(
1803 responseLibrary_->getResponse<panzer::Traits::Tangent>(responseName));
1804 resp->setVector(vec);
1805 is_active = true;
1806 }
1807 }
1808
1809 if (is_active) {
1810 for (std::size_t k=0; k<parameters_[j]->scalar_value.size(); k++) {
1811 std::string name = "PARAMETER_SENSITIVIES: "+(*parameters_[j]->names)[k];
1812 activeParameterNames.push_back(name);
1813 totalParameterCount++;
1814 }
1815 activeParameters.push_back(j);
1816 }
1817 }
1818
1819 // setup all the assembly in arguments
1821 setupAssemblyInArgs(inArgs,ae_inargs);
1822
1823 // add active parameter names to assembly in-args
1824 RCP<panzer::GlobalEvaluationData> ged_activeParameters =
1825 rcp(new panzer::ParameterList_GlobalEvaluationData(activeParameterNames));
1826 ae_inargs.addGlobalEvaluationData("PARAMETER_NAMES",ged_activeParameters);
1827
1828 // Initialize Fad components of all active parameters
1829 int paramIndex = 0;
1830 for (std::size_t ap=0; ap<activeParameters.size(); ++ap) {
1831 const int j = activeParameters[ap];
1832 for (unsigned int k=0; k < parameters_[j]->scalar_value.size(); k++) {
1833 panzer::Traits::FadType p(totalParameterCount, parameters_[j]->scalar_value[k].baseValue);
1834 p.fastAccessDx(paramIndex) = 1.0;
1835 parameters_[j]->scalar_value[k].family->template setValue<panzer::Traits::Tangent>(p);
1836 paramIndex++;
1837 }
1838 }
1839
1840 // make sure that the total parameter count and the total parameter index match up
1841 TEUCHOS_ASSERT(paramIndex==totalParameterCount);
1842
1843 // evaluate response tangent
1844 if(totalParameterCount>0) {
1845 responseLibrary_->addResponsesToInArgs<Traits::Tangent>(ae_inargs);
1846 responseLibrary_->evaluate<Traits::Tangent>(ae_inargs);
1847 }
1848}
1849
1850template <typename Scalar>
1851void
1853evalModelImpl_basic_dgdp_distro(const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
1854 const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
1855{
1856 typedef Thyra::ModelEvaluatorBase MEB;
1857
1858 // optional sanity check
1859 TEUCHOS_ASSERT(required_basic_dgdp_distro(outArgs));
1860
1861 // loop over parameters, and then build a dfdp_rl only if they are distributed
1862 // and the user has provided the UGI. Note that this may be overly expensive if they
1863 // don't actually want those sensitivites because memory will be allocated unneccesarily.
1864 // It would be good to do this "just in time", but for now this is sufficient.
1865 for(std::size_t p=0;p<parameters_.size();p++) {
1866
1867 // parameter is not distributed, a different path is
1868 // taken for those to compute dfdp
1869 if(!parameters_[p]->is_distributed)
1870 continue;
1871
1872 ResponseLibrary<Traits> & rLibrary = *parameters_[p]->dgdp_rl;
1873
1874 for(std::size_t r=0;r<responses_.size();r++) {
1875 // have derivatives been requested?
1876 MEB::Derivative<Scalar> deriv = outArgs.get_DgDp(r,p);
1877 if(deriv.isEmpty())
1878 continue;
1879
1880 Teuchos::RCP<Thyra::MultiVectorBase<Scalar> > vec = deriv.getMultiVector();
1881
1882 if(vec!=Teuchos::null) {
1883
1884 // get the response and tell it to fill the derivative vector
1885 std::string responseName = responses_[r]->name;
1886 Teuchos::RCP<panzer::ResponseMESupportBase<panzer::Traits::Jacobian> > resp
1887 = Teuchos::rcp_dynamic_cast<panzer::ResponseMESupportBase<panzer::Traits::Jacobian> >(
1888 rLibrary.getResponse<panzer::Traits::Jacobian>(responseName));
1889
1890 resp->setDerivative(vec);
1891 }
1892 }
1893
1894 // setup all the assembly in arguments (this is parameters and x/x_dot).
1895 // make sure the correct seeding is performed
1897 setupAssemblyInArgs(inArgs,ae_inargs);
1898
1899 ae_inargs.first_sensitivities_name = (*parameters_[p]->names)[0]; // distributed parameters have one name!
1900 ae_inargs.gather_seeds.push_back(1.0); // this assumes that gather point is always the zero index of
1901 // gather seeds
1902
1903 // evaluate responses
1904 rLibrary.addResponsesToInArgs<Traits::Jacobian>(ae_inargs);
1905 rLibrary.evaluate<Traits::Jacobian>(ae_inargs);
1906 }
1907}
1908
1909template <typename Scalar>
1910void
1912evalModelImpl_basic_dfdp_scalar(const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
1913 const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
1914{
1915 using Teuchos::RCP;
1916 using Teuchos::rcp_dynamic_cast;
1917
1918 typedef Thyra::ModelEvaluatorBase MEB;
1919
1920 TEUCHOS_ASSERT(required_basic_dfdp_scalar(outArgs));
1921
1922 // setup all the assembly in arguments (this is parameters and
1923 // x/x_dot). At this point with the exception of the one time dirichlet
1924 // beta that is all thats neccessary.
1926 setupAssemblyInArgs(inArgs,ae_inargs);
1927
1928 // First: Fill the output vectors from the input arguments structure. Put them
1929 // in the global evaluation data container so they are correctly communicated.
1931
1932 std::vector<std::string> activeParameters;
1933
1934 int totalParameterCount = 0;
1935 for(std::size_t i=0; i < parameters_.size(); i++) {
1936 // skip non-scalar parameters
1937 if(parameters_[i]->is_distributed)
1938 continue;
1939
1940 // have derivatives been requested?
1941 MEB::Derivative<Scalar> deriv = outArgs.get_DfDp(i);
1942 if(deriv.isEmpty())
1943 continue;
1944
1945 // grab multivector, make sure its the right dimension
1946 Teuchos::RCP<Thyra::MultiVectorBase<Scalar> > mVec = deriv.getMultiVector();
1947 TEUCHOS_ASSERT(mVec->domain()->dim()==Teuchos::as<int>(parameters_[i]->scalar_value.size()));
1948
1949 for (std::size_t j=0; j < parameters_[i]->scalar_value.size(); j++) {
1950
1951 // build containers for each vector
1952 RCP<LOCPair_GlobalEvaluationData> loc_pair
1953 = Teuchos::rcp(new LOCPair_GlobalEvaluationData(lof_,LinearObjContainer::F));
1954 RCP<LinearObjContainer> globalContainer = loc_pair->getGlobalLOC();
1955
1956 // stuff target vector into global container
1957 RCP<Thyra::VectorBase<Scalar> > vec = mVec->col(j);
1958 RCP<panzer::ThyraObjContainer<Scalar> > thGlobalContainer =
1959 Teuchos::rcp_dynamic_cast<panzer::ThyraObjContainer<Scalar> >(globalContainer);
1960 thGlobalContainer->set_f_th(vec);
1961
1962 // add container into in args object
1963 std::string name = "PARAMETER_SENSITIVIES: "+(*parameters_[i]->names)[j];
1964 ae_inargs.addGlobalEvaluationData(name,loc_pair->getGhostedLOC());
1965 ae_inargs.addGlobalEvaluationData(name+"_pair",loc_pair);
1966
1967 activeParameters.push_back(name);
1968 totalParameterCount++;
1969 }
1970 }
1971
1972 // Second: For all parameters that require derivative sensitivities, put in a name
1973 // so that the scatter can realize which sensitivity vectors it needs to fill
1975
1976 RCP<GlobalEvaluationData> ged_activeParameters
1977 = Teuchos::rcp(new ParameterList_GlobalEvaluationData(activeParameters));
1978 ae_inargs.addGlobalEvaluationData("PARAMETER_NAMES",ged_activeParameters);
1979
1980 // Third: Now seed all the parameters in the parameter vector so that derivatives
1981 // can be properly computed.
1983
1984 int paramIndex = 0;
1985 for(std::size_t i=0; i < parameters_.size(); i++) {
1986 // skip non-scalar parameters
1987 if(parameters_[i]->is_distributed)
1988 continue;
1989
1990 // don't modify the parameter if its not needed
1991 MEB::Derivative<Scalar> deriv = outArgs.get_DfDp(i);
1992 if(deriv.isEmpty()) {
1993 // reinitialize values that should not have sensitivities computed (this is a precaution)
1994 for (unsigned int j=0; j < parameters_[i]->scalar_value.size(); j++) {
1995 Traits::FadType p = Traits::FadType(totalParameterCount,
1996 parameters_[i]->scalar_value[j].baseValue);
1997 parameters_[i]->scalar_value[j].family->template setValue<panzer::Traits::Tangent>(p);
1998 }
1999 continue;
2000 }
2001 else {
2002 // loop over each parameter in the vector, initializing the AD type
2003 for (unsigned int j=0; j < parameters_[i]->scalar_value.size(); j++) {
2004 Traits::FadType p = Traits::FadType(totalParameterCount,
2005 parameters_[i]->scalar_value[j].baseValue);
2006 p.fastAccessDx(paramIndex) = 1.0;
2007 parameters_[i]->scalar_value[j].family->template setValue<panzer::Traits::Tangent>(p);
2008 paramIndex++;
2009 }
2010 }
2011 }
2012
2013 // make sure that the total parameter count and the total parameter index match up
2014 TEUCHOS_ASSERT(paramIndex==totalParameterCount);
2015
2016 // Fourth: Actually evaluate the residual's sensitivity to the parameters
2018
2019 if(totalParameterCount>0) {
2020 PANZER_FUNC_TIME_MONITOR("panzer::ModelEvaluator::evalModel(df/dp)");
2021 ae_tm_.getAsObject<panzer::Traits::Tangent>()->evaluate(ae_inargs);
2022 }
2023}
2024
2025template <typename Scalar>
2026void
2028evalModelImpl_basic_dfdp_scalar_fd(const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
2029 const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
2030{
2031 PANZER_FUNC_TIME_MONITOR("panzer::ModelEvaluator::evalModel(df/dp)");
2032
2033 using Teuchos::RCP;
2034 using Teuchos::rcp_dynamic_cast;
2035
2036 typedef Thyra::ModelEvaluatorBase MEB;
2037
2038 TEUCHOS_ASSERT(required_basic_dfdp_scalar(outArgs));
2039
2040 // First evaluate the model without df/dp for the base point
2041 // Maybe it would be better to set all outArgs and then remove the df/dp ones,
2042 // but I couldn't get that to work.
2043 MEB::OutArgs<Scalar> outArgs_base = this->createOutArgs();
2044 if (outArgs.get_f() == Teuchos::null)
2045 outArgs_base.set_f(Thyra::createMember(this->get_f_space()));
2046 else
2047 outArgs_base.set_f(outArgs.get_f());
2048 outArgs_base.set_W_op(outArgs.get_W_op());
2049 this->evalModel(inArgs, outArgs_base);
2050 RCP<const Thyra::VectorBase<Scalar> > f = outArgs_base.get_f();
2051 RCP<const Thyra::VectorBase<Scalar> > x = inArgs.get_x();
2052 RCP<const Thyra::VectorBase<Scalar> > x_dot;
2053 if (inArgs.supports(MEB::IN_ARG_x_dot))
2054 x_dot = inArgs.get_x_dot();
2055
2056 // Create in/out args for FD calculation
2057 RCP<Thyra::VectorBase<Scalar> > fd = Thyra::createMember(this->get_f_space());
2058 MEB::OutArgs<Scalar> outArgs_fd = this->createOutArgs();
2059 outArgs_fd.set_f(fd);
2060
2061 RCP<Thyra::VectorBase<Scalar> > xd = Thyra::createMember(this->get_x_space());
2062 RCP<Thyra::VectorBase<Scalar> > xd_dot;
2063 if (x_dot != Teuchos::null)
2064 xd_dot = Thyra::createMember(this->get_x_space());
2065 MEB::InArgs<Scalar> inArgs_fd = this->createInArgs();
2066 inArgs_fd.setArgs(inArgs); // This sets all inArgs that we don't override below
2067 inArgs_fd.set_x(xd);
2068 if (x_dot != Teuchos::null)
2069 inArgs_fd.set_x_dot(xd_dot);
2070
2071 const double h = fd_perturb_size_;
2072 for(std::size_t i=0; i < parameters_.size(); i++) {
2073
2074 // skip non-scalar parameters
2075 if(parameters_[i]->is_distributed)
2076 continue;
2077
2078 // have derivatives been requested?
2079 MEB::Derivative<Scalar> deriv = outArgs.get_DfDp(i);
2080 if(deriv.isEmpty())
2081 continue;
2082
2083 // grab multivector, make sure its the right dimension
2084 RCP<Thyra::MultiVectorBase<Scalar> > dfdp = deriv.getMultiVector();
2085 TEUCHOS_ASSERT(dfdp->domain()->dim()==Teuchos::as<int>(parameters_[i]->scalar_value.size()));
2086
2087 // Get parameter vector and tangent vectors
2088 RCP<const Thyra::VectorBase<Scalar> > p = inArgs.get_p(i);
2089 RCP<const Thyra::VectorBase<Scalar> > dx_v = inArgs.get_p(i+parameters_.size());
2090 RCP<const Thyra::MultiVectorBase<Scalar> > dx =
2091 rcp_dynamic_cast<const Thyra::DefaultMultiVectorProductVector<Scalar> >(dx_v,true)->getMultiVector();
2092 RCP<const Thyra::VectorBase<Scalar> > dx_dot_v;
2093 RCP<const Thyra::MultiVectorBase<Scalar> > dx_dot;
2094 if (x_dot != Teuchos::null) {
2095 dx_dot_v =inArgs.get_p(i+parameters_.size()+tangent_space_.size());
2096 dx_dot =
2097 rcp_dynamic_cast<const Thyra::DefaultMultiVectorProductVector<Scalar> >(dx_dot_v,true)->getMultiVector();
2098 }
2099
2100 // Create perturbed parameter vector
2101 RCP<Thyra::VectorBase<Scalar> > pd = Thyra::createMember(this->get_p_space(i));
2102 inArgs_fd.set_p(i,pd);
2103
2104 for (std::size_t j=0; j < parameters_[i]->scalar_value.size(); j++) {
2105
2106 // Perturb parameter vector
2107 Thyra::copy(*p, pd.ptr());
2108 Thyra::set_ele(j, Thyra::get_ele(*p,j)+h, pd.ptr());
2109
2110 // Perturb state vectors using tangents
2111 Thyra::V_VpStV(xd.ptr(), *x, h, *(dx)->col(j));
2112 if (x_dot != Teuchos::null)
2113 Thyra::V_VpStV(xd_dot.ptr(), *x_dot, h, *(dx_dot)->col(j));
2114
2115 // Evaluate perturbed residual
2116 Thyra::assign(fd.ptr(), 0.0);
2117 this->evalModel(inArgs_fd, outArgs_fd);
2118
2119 // FD calculation
2120 Thyra::V_StVpStV(dfdp->col(j).ptr(), 1.0/h, *fd, -1.0/h, *f);
2121
2122 // Reset parameter back to un-perturbed value
2123 parameters_[i]->scalar_value[j].family->setRealValueForAllTypes(Thyra::get_ele(*p,j));
2124
2125 }
2126 }
2127}
2128
2129template <typename Scalar>
2130void
2132evalModelImpl_basic_dfdp_distro(const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
2133 const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
2134{
2135 using Teuchos::RCP;
2136 using Teuchos::rcp_dynamic_cast;
2137 using Teuchos::null;
2138
2139 typedef Thyra::ModelEvaluatorBase MEB;
2140
2141 TEUCHOS_ASSERT(required_basic_dfdp_distro(outArgs));
2142
2143 // loop over parameters, and then build a dfdp_rl only if they are distributed
2144 // and the user has provided the UGI. Note that this may be overly expensive if they
2145 // don't actually want those sensitivites because memory will be allocated unneccesarily.
2146 // It would be good to do this "just in time", but for now this is sufficient.
2147 for(std::size_t p=0;p<parameters_.size();p++) {
2148
2149 // parameter is not distributed, a different path is
2150 // taken for those to compute dfdp
2151 if(!parameters_[p]->is_distributed)
2152 continue;
2153
2154 // parameter is distributed but has no global indexer.
2155 // thus the user doesn't want sensitivities!
2156 if(parameters_[p]->dfdp_rl==null)
2157 continue;
2158
2159 // have derivatives been requested?
2160 MEB::Derivative<Scalar> deriv = outArgs.get_DfDp(p);
2161 if(deriv.isEmpty())
2162 continue;
2163
2164 ResponseLibrary<Traits> & rLibrary = *parameters_[p]->dfdp_rl;
2165
2166 // get the response and tell it to fill the derivative operator
2167 RCP<Response_Residual<Traits::Jacobian> > response_jacobian =
2168 rcp_dynamic_cast<Response_Residual<Traits::Jacobian> >(rLibrary.getResponse<Traits::Jacobian>("RESIDUAL"));
2169 response_jacobian->setJacobian(deriv.getLinearOp());
2170
2171 // setup all the assembly in arguments (this is parameters and x/x_dot).
2172 // make sure the correct seeding is performed
2174 setupAssemblyInArgs(inArgs,ae_inargs);
2175
2176 ae_inargs.first_sensitivities_name = (*parameters_[p]->names)[0]; // distributed parameters have one name!
2177 ae_inargs.gather_seeds.push_back(1.0); // this assumes that gather point is always the zero index of
2178 // gather seeds
2179 rLibrary.addResponsesToInArgs<Traits::Jacobian>(ae_inargs);
2180
2181 rLibrary.evaluate<Traits::Jacobian>(ae_inargs);
2182 }
2183}
2184
2185template <typename Scalar>
2187required_basic_g(const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
2188{
2189 // determine if any of the outArgs are not null!
2190 bool activeGArgs = false;
2191 for(int i=0;i<outArgs.Ng();i++)
2192 activeGArgs |= (outArgs.get_g(i)!=Teuchos::null);
2193
2194 return activeGArgs | required_basic_dgdx(outArgs);
2195}
2196
2197template <typename Scalar>
2199required_basic_dgdx(const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
2200{
2201 typedef Thyra::ModelEvaluatorBase MEB;
2202
2203 // determine if any of the outArgs are not null!
2204 bool activeGArgs = false;
2205 for(int i=0;i<outArgs.Ng();i++) {
2206 // no derivatives are supported
2207 if(outArgs.supports(MEB::OUT_ARG_DgDx,i).none())
2208 continue;
2209
2210 // this is basically a redundant computation
2211 activeGArgs |= (!outArgs.get_DgDx(i).isEmpty());
2212 }
2213
2214 return activeGArgs;
2215}
2216
2217template <typename Scalar>
2219required_basic_dgdp_scalar(const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
2220{
2221 typedef Thyra::ModelEvaluatorBase MEB;
2222
2223 // determine if any of the outArgs are not null!
2224 bool activeGArgs = false;
2225 for(int i=0;i<outArgs.Ng();i++) {
2226 for(int p=0;p<Teuchos::as<int>(parameters_.size());p++) {
2227
2228 // only look at scalar parameters
2229 if(parameters_[p]->is_distributed)
2230 continue;
2231
2232 // no derivatives are supported
2233 if(outArgs.supports(MEB::OUT_ARG_DgDp,i,p).none())
2234 continue;
2235
2236 activeGArgs |= (!outArgs.get_DgDp(i,p).isEmpty());
2237 }
2238 }
2239
2240 return activeGArgs;
2241}
2242
2243template <typename Scalar>
2245required_basic_dgdp_distro(const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
2246{
2247 typedef Thyra::ModelEvaluatorBase MEB;
2248
2249 // determine if any of the outArgs are not null!
2250 bool activeGArgs = false;
2251 for(int i=0;i<outArgs.Ng();i++) {
2252 for(int p=0;p<Teuchos::as<int>(parameters_.size());p++) {
2253
2254 // only look at distributed parameters
2255 if(!parameters_[p]->is_distributed)
2256 continue;
2257
2258 // no derivatives are supported
2259 if(outArgs.supports(MEB::OUT_ARG_DgDp,i,p).none())
2260 continue;
2261
2262 activeGArgs |= (!outArgs.get_DgDp(i,p).isEmpty());
2263 }
2264 }
2265
2266 return activeGArgs;
2267}
2268
2269template <typename Scalar>
2271required_basic_dfdp_scalar(const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
2272{
2273 typedef Thyra::ModelEvaluatorBase MEB;
2274
2275 // determine if any of the outArgs are not null!
2276 bool activeFPArgs = false;
2277 for(int i=0;i<Teuchos::as<int>(parameters_.size());i++) {
2278
2279 // this is for scalar parameters only
2280 if(parameters_[i]->is_distributed)
2281 continue;
2282
2283 // no derivatives are supported
2284 if(outArgs.supports(MEB::OUT_ARG_DfDp,i).none())
2285 continue;
2286
2287 // this is basically a redundant computation
2288 activeFPArgs |= (!outArgs.get_DfDp(i).isEmpty());
2289 }
2290
2291 return activeFPArgs;
2292}
2293
2294template <typename Scalar>
2296required_basic_dfdp_distro(const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
2297{
2298 typedef Thyra::ModelEvaluatorBase MEB;
2299
2300 // determine if any of the outArgs are not null!
2301 bool activeFPArgs = false;
2302 for(int i=0;i<Teuchos::as<int>(parameters_.size());i++) {
2303
2304 // this is for scalar parameters only
2305 if(!parameters_[i]->is_distributed)
2306 continue;
2307
2308 // no derivatives are supported
2309 if(outArgs.supports(MEB::OUT_ARG_DfDp,i).none())
2310 continue;
2311
2312 // this is basically a redundant computation
2313 activeFPArgs |= (!outArgs.get_DfDp(i).isEmpty());
2314 }
2315
2316 return activeFPArgs;
2317}
2318
2319template <typename Scalar>
2322 const Teuchos::RCP<panzer::WorksetContainer> & wc,
2323 const std::vector<Teuchos::RCP<panzer::PhysicsBlock> >& physicsBlocks,
2324 const std::vector<panzer::BC> & bcs,
2325 const panzer::EquationSetFactory & eqset_factory,
2326 const panzer::BCStrategyFactory& bc_factory,
2328 const Teuchos::ParameterList& closure_models,
2329 const Teuchos::ParameterList& user_data,
2330 const bool write_graphviz_file,
2331 const std::string& graphviz_file_prefix)
2332{
2333 using Teuchos::RCP;
2334 using Teuchos::rcp;
2335 using Teuchos::null;
2336
2337 // loop over parameters, and then build a dfdp_rl only if they are distributed
2338 // and the user has provided the UGI. Note that this may be overly expensive if they
2339 // don't actually want those sensitivites because memory will be allocated unneccesarily.
2340 // It would be good to do this "just in time", but for now this is sufficient.
2341 for(std::size_t p=0;p<parameters_.size();p++) {
2342 // parameter is not distributed, a different path is
2343 // taken for those to compute dfdp
2344 if(!parameters_[p]->is_distributed)
2345 continue;
2346
2347 // parameter is distributed but has no global indexer.
2348 // thus the user doesn't want sensitivities!
2349 if(parameters_[p]->global_indexer==null)
2350 continue;
2351
2352 // build the linear object factory that has the correct sizing for
2353 // the sensitivity matrix (parameter sized domain, residual sized range)
2354 RCP<const LinearObjFactory<Traits> > param_lof = cloneWithNewDomain(*lof_,
2355 parameters_[p]->global_indexer);
2356
2357 // the user wants global sensitivities, hooray! Build and setup the response library
2358 RCP<ResponseLibrary<Traits> > rLibrary
2359 = Teuchos::rcp(new ResponseLibrary<Traits>(wc,lof_->getRangeGlobalIndexer(),
2360 param_lof,true));
2361 rLibrary->buildResidualResponseEvaluators(physicsBlocks,eqset_factory,bcs,bc_factory,
2362 cm_factory,closure_models,user_data,
2363 write_graphviz_file,graphviz_file_prefix);
2364
2365 // make sure parameter response library is correct
2366 parameters_[p]->dfdp_rl = rLibrary;
2367 }
2368}
2369
2370template <typename Scalar>
2373 const Teuchos::RCP<panzer::WorksetContainer> & wc,
2374 const std::vector<Teuchos::RCP<panzer::PhysicsBlock> >& physicsBlocks,
2375 const std::vector<panzer::BC>& /* bcs */,
2376 const panzer::EquationSetFactory & eqset_factory,
2377 const panzer::BCStrategyFactory& /* bc_factory */,
2379 const Teuchos::ParameterList& closure_models,
2380 const Teuchos::ParameterList& user_data,
2381 const bool write_graphviz_file,
2382 const std::string& graphviz_file_prefix)
2383{
2384 using Teuchos::RCP;
2385 using Teuchos::rcp;
2386 using Teuchos::null;
2387
2388 // loop over parameters, and then build a dfdp_rl only if they are distributed
2389 // and the user has provided the UGI. Note that this may be overly expensive if they
2390 // don't actually want those sensitivites because memory will be allocated unneccesarily.
2391 // It would be good to do this "just in time", but for now this is sufficient.
2392 for(std::size_t p=0;p<parameters_.size();p++) {
2393 // parameter is not distributed, a different path is
2394 // taken for those to compute dfdp
2395 if(!parameters_[p]->is_distributed)
2396 continue;
2397
2398 // parameter is distributed but has no global indexer.
2399 // thus the user doesn't want sensitivities!
2400 if(parameters_[p]->global_indexer==null)
2401 continue;
2402
2403 // extract the linear object factory that has the correct sizing for
2404 // the sensitivity vector
2405 RCP<const LinearObjFactory<Traits> > param_lof = parameters_[p]->dfdp_rl->getLinearObjFactory();
2406 RCP<const GlobalIndexer > param_ugi = parameters_[p]->global_indexer;
2407
2408 // the user wants global sensitivities, hooray! Build and setup the response library
2409 RCP<ResponseLibrary<Traits> > rLibrary
2410 = Teuchos::rcp(new ResponseLibrary<Traits>(wc,lof_->getRangeGlobalIndexer(), lof_));
2411
2412
2413 // build evaluators for all flexible responses
2414 for(std::size_t r=0;r<responses_.size();r++) {
2415 // only responses with a builder are non null!
2416 if(responses_[r]->builder==Teuchos::null)
2417 continue;
2418
2419 // set the current derivative information in the builder
2420 // responses_[r]->builder->setDerivativeInformationBase(param_lof,param_ugi);
2421 responses_[r]->builder->setDerivativeInformation(param_lof);
2422
2423 // add the response
2424 rLibrary->addResponse(responses_[r]->name,
2425 responses_[r]->wkst_desc,
2426 *responses_[r]->builder);
2427 }
2428
2429 rLibrary->buildResponseEvaluators(physicsBlocks,eqset_factory,
2430 cm_factory,closure_models,user_data,
2431 write_graphviz_file,graphviz_file_prefix);
2432
2433 // make sure parameter response library is correct
2434 parameters_[p]->dgdp_rl = rLibrary;
2435 }
2436}
2437
2438template <typename Scalar>
2440setOneTimeDirichletBeta(const Scalar & beta) const
2441{
2442 oneTimeDirichletBeta_on_ = true;
2443 oneTimeDirichletBeta_ = beta;
2444}
2445
2446template <typename Scalar>
2447Teuchos::RCP<typename panzer::ModelEvaluator<Scalar>::ParameterObject>
2449createScalarParameter(const Teuchos::Array<std::string> & in_names,
2450 const Teuchos::Array<Scalar> & in_values) const
2451{
2452 using Teuchos::RCP;
2453 using Teuchos::rcp;
2454 using Teuchos::rcp_dynamic_cast;
2455 using Teuchos::ptrFromRef;
2456
2457 TEUCHOS_ASSERT(in_names.size()==in_values.size());
2458
2459 // Check that the parameters are valid (i.e., they already exist in the parameter library)
2460 // std::size_t np = in_names.size();
2461 // for(std::size_t i=0;i<np;i++)
2462 // TEUCHOS_TEST_FOR_EXCEPTION(!global_data_->pl->isParameter(in_names[i]),
2463 // std::logic_error,
2464 // "Parameter \"" << in_names[i] << "\" does not exist in parameter library!");
2465
2466 RCP<ParameterObject> paramObj = rcp(new ParameterObject);
2467
2468 paramObj->names = rcp(new Teuchos::Array<std::string>(in_names));
2469 paramObj->is_distributed = false;
2470
2471 // register all the scalar parameters, setting initial
2472 for(int i=0;i<in_names.size();i++)
2473 registerScalarParameter(in_names[i],*global_data_->pl,in_values[i]);
2474
2475 paramObj->scalar_value = panzer::ParamVec();
2476 global_data_->pl->fillVector<panzer::Traits::Residual>(*paramObj->names, paramObj->scalar_value);
2477
2478 // build initial condition vector
2479 paramObj->space =
2480 Thyra::locallyReplicatedDefaultSpmdVectorSpace<Scalar>(
2481 rcp(new Teuchos::MpiComm<long int>(lof_->getComm().getRawMpiComm())),paramObj->names->size());
2482
2483 // fill vector with parameter values
2484 Teuchos::ArrayRCP<Scalar> data;
2485 RCP<Thyra::VectorBase<Scalar> > initial_value = Thyra::createMember(paramObj->space);
2486 RCP<Thyra::SpmdVectorBase<Scalar> > vec = rcp_dynamic_cast<Thyra::SpmdVectorBase<Scalar> >(initial_value);
2487 vec->getNonconstLocalData(ptrFromRef(data));
2488 for (unsigned int i=0; i < paramObj->scalar_value.size(); i++)
2489 data[i] = in_values[i];
2490
2491 paramObj->initial_value = initial_value;
2492
2493 return paramObj;
2494}
2495
2496template <typename Scalar>
2497Teuchos::RCP<typename panzer::ModelEvaluator<Scalar>::ParameterObject>
2499createDistributedParameter(const std::string & key,
2500 const Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > & vs,
2501 const Teuchos::RCP<const Thyra::VectorBase<Scalar> > & initial,
2502 const Teuchos::RCP<const GlobalIndexer> & ugi) const
2503{
2504 using Teuchos::RCP;
2505 using Teuchos::rcp;
2506
2507 RCP<ParameterObject> paramObj = rcp(new ParameterObject);
2508
2509 paramObj->is_distributed = true;
2510 paramObj->names = rcp(new Teuchos::Array<std::string>());
2511 paramObj->names->push_back(key);
2512 paramObj->space = vs;
2513 paramObj->initial_value = initial;
2514
2515 paramObj->global_indexer = ugi;
2516
2517 return paramObj;
2518}
2519
2520template <typename Scalar>
2521void
2523setParameters(const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs) const
2524{
2525 for(std::size_t i=0; i < parameters_.size(); i++) {
2526
2527 // skip non-scalar parameters (for now)
2528 if(parameters_[i]->is_distributed)
2529 continue;
2530
2531 // set parameter values for given parameter vector for all evaluation types
2532 Teuchos::RCP<const Thyra::VectorBase<Scalar> > p = inArgs.get_p(i);
2533 if (p != Teuchos::null) {
2534 for (unsigned int j=0; j < parameters_[i]->scalar_value.size(); j++) {
2535 parameters_[i]->scalar_value[j].family->setRealValueForAllTypes(Thyra::get_ele(*p,j));
2536 }
2537 }
2538
2539 }
2540}
2541
2542template <typename Scalar>
2543void
2545resetParameters() const
2546{
2547 for(std::size_t i=0; i < parameters_.size(); i++) {
2548
2549 // skip non-scalar parameters (for now)
2550 if(parameters_[i]->is_distributed)
2551 continue;
2552
2553 // Reset each parameter back to its nominal
2554 for (unsigned int j=0; j < parameters_[i]->scalar_value.size(); j++) {
2555 parameters_[i]->scalar_value[j].family->setRealValueForAllTypes(Thyra::get_ele(*(parameters_[i]->initial_value),j));
2556 }
2557
2558 }
2559}
2560
2561#endif // __Panzer_ModelEvaluator_impl_hpp__
PHX::MDField< ScalarT, panzer::Cell, panzer::IP > result
A field that will be used to build up the result of the integral we're performing.
void addGlobalEvaluationData(const std::string &key, const Teuchos::RCP< GlobalEvaluationData > &ged)
Teuchos::RCP< panzer::LinearObjContainer > ghostedContainer_
Teuchos::RCP< panzer::LinearObjContainer > container_
virtual void evalModelImpl_basic_dfdp_scalar_fd(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const
void setOneTimeDirichletBeta(const Scalar &beta) const
virtual void evalModelImpl_basic_dgdp_scalar(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const
void setupAssemblyInArgs(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, panzer::AssemblyEngineInArgs &ae_inargs) const
void setupModel(const Teuchos::RCP< panzer::WorksetContainer > &wc, const std::vector< Teuchos::RCP< panzer::PhysicsBlock > > &physicsBlocks, const std::vector< panzer::BC > &bcs, const panzer::EquationSetFactory &eqset_factory, const panzer::BCStrategyFactory &bc_factory, const panzer::ClosureModelFactory_TemplateManager< panzer::Traits > &volume_cm_factory, const panzer::ClosureModelFactory_TemplateManager< panzer::Traits > &bc_cm_factory, const Teuchos::ParameterList &closure_models, const Teuchos::ParameterList &user_data, bool writeGraph=false, const std::string &graphPrefix="", const Teuchos::ParameterList &me_params=Teuchos::ParameterList())
Thyra::ModelEvaluatorBase::OutArgs< Scalar > createOutArgsImpl() const override
bool required_basic_g(const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const
Does this set of out args require a simple response?
void evalModel_D2fDp2(int pIndex, const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Teuchos::RCP< const Thyra::VectorBase< Scalar > > &delta_x, const Teuchos::RCP< Thyra::LinearOpBase< Scalar > > &D2fDp2) const
Teuchos::RCP< Thyra::LinearOpBase< Scalar > > create_DfDp_op(int i) const override
virtual void evalModelImpl_basic_dgdx(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > f_space_
void buildDistroParamDfDp_RL(const Teuchos::RCP< panzer::WorksetContainer > &wc, const std::vector< Teuchos::RCP< panzer::PhysicsBlock > > &physicsBlocks, const std::vector< panzer::BC > &bcs, const panzer::EquationSetFactory &eqset_factory, const panzer::BCStrategyFactory &bc_factory, const panzer::ClosureModelFactory_TemplateManager< panzer::Traits > &cm_factory, const Teuchos::ParameterList &closure_models, const Teuchos::ParameterList &user_data, const bool write_graphviz_file=false, const std::string &graphviz_file_prefix="")
Teuchos::RCP< const Thyra::LinearOpWithSolveFactoryBase< Scalar > > get_W_factory() const override
void setParameters(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs) const
Teuchos::RCP< ParameterObject > createDistributedParameter(const std::string &key, const Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > &vs, const Teuchos::RCP< const Thyra::VectorBase< Scalar > > &initial, const Teuchos::RCP< const GlobalIndexer > &ugi) const
panzer::AssemblyEngine_TemplateManager< panzer::Traits > ae_tm_
virtual void evalModelImpl_basic_dfdp_scalar(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const
bool required_basic_dfdp_scalar(const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const
Are derivatives of the residual with respect to the scalar parameters in the out args?...
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > x_space_
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_p_space(int i) const override
virtual void evalModelImpl_basic_g(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const
Construct a simple response dicatated by this set of out args.
void evalModel_D2fDx2(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Teuchos::RCP< const Thyra::VectorBase< Scalar > > &delta_x, const Teuchos::RCP< Thyra::LinearOpBase< Scalar > > &D2fDx2) const
bool required_basic_dgdx(const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const
Are their required responses in the out args? DgDx.
bool required_basic_dgdp_distro(const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const
Are their required responses in the out args? DgDp.
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() const override
virtual void evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const override
void evalModel_D2fDpDx(int pIndex, const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Teuchos::RCP< const Thyra::VectorBase< Scalar > > &delta_x, const Teuchos::RCP< Thyra::LinearOpBase< Scalar > > &D2fDpDx) const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_f_space() const override
void buildVolumeFieldManagers(const bool value)
virtual void evalModelImpl_basic_dgdp_distro(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const
void evalModel_D2gDxDp(int rIndex, int pIndex, const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Teuchos::RCP< const Thyra::VectorBase< Scalar > > &delta_p, const Teuchos::RCP< Thyra::VectorBase< Scalar > > &D2gDxDp) const
int addDistributedParameter(const std::string &name, const Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > &vs, const Teuchos::RCP< GlobalEvaluationData > &ged, const Teuchos::RCP< const Thyra::VectorBase< Scalar > > &initial, const Teuchos::RCP< const GlobalIndexer > &ugi=Teuchos::null)
bool required_basic_dgdp_scalar(const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const
Are their required responses in the out args? DgDp.
int addParameter(const std::string &name, const Scalar &initial)
Thyra::ModelEvaluatorBase::InArgs< Scalar > getNominalValues() const override
const std::string & get_g_name(int i) const
void buildDistroParamDgDp_RL(const Teuchos::RCP< panzer::WorksetContainer > &wc, const std::vector< Teuchos::RCP< panzer::PhysicsBlock > > &physicsBlocks, const std::vector< panzer::BC > &bcs, const panzer::EquationSetFactory &eqset_factory, const panzer::BCStrategyFactory &bc_factory, const panzer::ClosureModelFactory_TemplateManager< panzer::Traits > &cm_factory, const Teuchos::ParameterList &closure_models, const Teuchos::ParameterList &user_data, const bool write_graphviz_file=false, const std::string &graphviz_file_prefix="")
Teuchos::RCP< const Teuchos::Array< std::string > > get_p_names(int i) const override
Teuchos::ArrayView< const std::string > get_g_names(int i) const override
void applyDirichletBCs(const Teuchos::RCP< Thyra::VectorBase< Scalar > > &x, const Teuchos::RCP< Thyra::VectorBase< Scalar > > &f) const
Thyra::ModelEvaluatorBase::InArgs< Scalar > createInArgs() const override
int addFlexibleResponse(const std::string &responseName, const std::vector< WorksetDescriptor > &wkst_desc, const Teuchos::RCP< ResponseMESupportBuilderBase > &builder)
void evalModel_D2gDpDx(int rIndex, int pIndex, const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Teuchos::RCP< const Thyra::VectorBase< Scalar > > &delta_x, const Teuchos::RCP< Thyra::VectorBase< Scalar > > &D2gDpDx) const
void buildBCFieldManagers(const bool value)
Teuchos::RCP< const panzer::LinearObjFactory< panzer::Traits > > lof_
virtual void evalModelImpl_basic_dfdp_distro(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const
void evalModel_D2fDxDp(int pIndex, const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Teuchos::RCP< const Thyra::VectorBase< Scalar > > &delta_p, const Teuchos::RCP< Thyra::LinearOpBase< Scalar > > &D2fDxDp) const
void evalModel_D2gDx2(int rIndex, const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Teuchos::RCP< const Thyra::VectorBase< Scalar > > &delta_x, const Teuchos::RCP< Thyra::VectorBase< Scalar > > &D2gDx2) const
virtual void evalModelImpl_basic(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const
Evaluate a simple model, meaning a residual and a jacobian, no fancy stochastic galerkin or multipoin...
Teuchos::RCP< Thyra::LinearOpBase< Scalar > > create_W_op() const override
Teuchos::RCP< panzer::ResponseLibrary< panzer::Traits > > responseLibrary_
void evalModel_D2gDp2(int rIndex, int pIndex, const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Teuchos::RCP< const Thyra::VectorBase< Scalar > > &delta_x, const Teuchos::RCP< Thyra::VectorBase< Scalar > > &D2gDp2) const
bool required_basic_dfdp_distro(const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const
Are derivatives of the residual with respect to the distributed parameters in the out args?...
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_g_space(int i) const override
void addNonParameterGlobalEvaluationData(const std::string &name, const Teuchos::RCP< GlobalEvaluationData > &ged)
void initializeNominalValues() const
Initialize the nominal values with good starting conditions.
Teuchos::RCP< ParameterObject > createScalarParameter(const Teuchos::Array< std::string > &names, const Teuchos::Array< Scalar > &in_values) const
Teuchos::RCP< ResponseBase > getResponse(const std::string &responseName) const
void evaluate(const panzer::AssemblyEngineInArgs &input_args)
void addResponsesToInArgs(panzer::AssemblyEngineInArgs &input_args) const
Teuchos::RCP< const LinearObjFactory< panzer::Traits > > cloneWithNewDomain(const LinearObjFactory< panzer::Traits > &lof, const Teuchos::RCP< const GlobalIndexer > &dUgi)
Clone a linear object factory, but using a different domain.
Sacado::ScalarParameterVector< panzer::EvaluationTraits > ParamVec
Kokkos::Compat::KokkosDeviceWrapperNode< PHX::Device > TpetraNodeType
void registerScalarParameter(const std::string name, panzer::ParamLib &pl, double realValue)
Interface for constructing a BCStrategy_TemplateManager.
Allocates and initializes an equation set template manager.
Teuchos::RCP< panzer::ResponseLibrary< panzer::Traits > > dfdp_rl
Teuchos::RCP< const GlobalIndexer > global_indexer
PANZER_FADTYPE FadType