Panzer Version of the Day
Loading...
Searching...
No Matches
Panzer_ModelEvaluator_Epetra.cpp
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#include "PanzerDiscFE_config.hpp"
50#include "Panzer_GlobalData.hpp"
54
55#include "Epetra_CrsMatrix.h"
56#include "Epetra_MpiComm.h"
57#include "Epetra_LocalMap.h"
58
59#include "Teuchos_ScalarTraits.hpp"
60#include "Teuchos_DefaultMpiComm.hpp"
61#include "Teuchos_TimeMonitor.hpp"
62
63#include "Thyra_EpetraOperatorWrapper.hpp"
64#include "Thyra_EpetraThyraWrappers.hpp"
65#include "Thyra_VectorStdOps.hpp"
66#include "Thyra_ProductVectorBase.hpp"
67#include "Thyra_SpmdVectorBase.hpp"
68#include "Thyra_DefaultProductVector.hpp"
69#include "Thyra_ProductVectorSpaceBase.hpp"
70
71#include <sstream>
72
73namespace {
74
76class EpetraLOF_EOpFactory : public Teuchos::AbstractFactory<Epetra_Operator> {
77 Teuchos::RCP<Epetra_CrsGraph> W_graph_;
78public:
79 EpetraLOF_EOpFactory(const panzer::BlockedEpetraLinearObjFactory<panzer::Traits,int> & lof)
80 : W_graph_(lof.getGraph(0,0)) {}
81
82 virtual Teuchos::RCP<Epetra_Operator> create() const
83 { return Teuchos::rcp(new Epetra_CrsMatrix(::Copy,*W_graph_)); }
84};
85
86}
87
88
90ModelEvaluator_Epetra(const Teuchos::RCP<panzer::FieldManagerBuilder>& fmb,
91 const Teuchos::RCP<panzer::ResponseLibrary<panzer::Traits> >& rLibrary,
92 const Teuchos::RCP<panzer::LinearObjFactory<panzer::Traits> > & lof,
93 const std::vector<Teuchos::RCP<Teuchos::Array<std::string> > >& p_names,
94 const std::vector<Teuchos::RCP<Teuchos::Array<double> > >& p_values,
95 const Teuchos::RCP<panzer::GlobalData>& global_data,
96 bool build_transient_support)
97 : t_init_(0.0)
98 , fmb_(fmb)
99 , responseLibrary_(rLibrary)
100 , p_names_(p_names)
101 , global_data_(global_data)
102 , build_transient_support_(build_transient_support)
103 , lof_(lof)
104 , oneTimeDirichletBeta_on_(false)
105 , oneTimeDirichletBeta_(0.0)
106{
107 using Teuchos::rcp;
108 using Teuchos::rcp_dynamic_cast;
109
111 ae_tm_.buildObjects(builder);
112
113 // Setup parameters
114 this->initializeParameterVector(p_names_,p_values,global_data->pl);
115
116 // try to determine the runtime linear object factory
117
118 Teuchos::RCP<panzer::BlockedEpetraLinearObjFactory<panzer::Traits,int> > ep_lof =
119 Teuchos::rcp_dynamic_cast<panzer::BlockedEpetraLinearObjFactory<panzer::Traits,int> >(lof);
120
121 // initialize maps, x_dot_init, x0, p_init, g_map, and linear operator factory
122 if(ep_lof!=Teuchos::null)
123 initializeEpetraObjs(*ep_lof);
124 else {
125 TEUCHOS_ASSERT(false); // bad news!
126 }
127}
128
130ModelEvaluator_Epetra(const Teuchos::RCP<panzer::FieldManagerBuilder>& fmb,
131 const Teuchos::RCP<panzer::ResponseLibrary<panzer::Traits> >& rLibrary,
133 const std::vector<Teuchos::RCP<Teuchos::Array<std::string> > >& p_names,
134 const std::vector<Teuchos::RCP<Teuchos::Array<double> > >& p_values,
135 const Teuchos::RCP<panzer::GlobalData>& global_data,
136 bool build_transient_support)
137 : t_init_(0.0)
138 , fmb_(fmb)
139 , responseLibrary_(rLibrary)
140 , p_names_(p_names)
141 , global_data_(global_data)
142 , build_transient_support_(build_transient_support)
143 , lof_(lof)
144 , oneTimeDirichletBeta_on_(false)
145 , oneTimeDirichletBeta_(0.0)
146{
147 using Teuchos::rcp;
148 using Teuchos::rcp_dynamic_cast;
149
151 ae_tm_.buildObjects(builder);
152
153 // Setup parameters
154 this->initializeParameterVector(p_names_,p_values,global_data->pl);
155
156 // initailize maps, x_dot_init, x0, p_init, g_map, and linear operator factory
158}
159
161{
162 using Teuchos::RCP;
163 using Teuchos::rcp;
164 using Teuchos::rcp_dynamic_cast;
165
166 TEUCHOS_TEST_FOR_EXCEPTION(responseLibrary_==Teuchos::null,std::logic_error,
167 "panzer::ModelEvaluator_Epetra::initializeEpetraObjs: The response library "
168 "was not correctly initialized before calling initializeEpetraObjs.");
169
170 map_x_ = lof.getMap(0);
171 x0_ = rcp(new Epetra_Vector(*map_x_));
172 x_dot_init_ = rcp(new Epetra_Vector(*map_x_));
173 x_dot_init_->PutScalar(0.0);
174
175 // setup parameters (initialize vector from parameter library)
176 for (int i=0; i < parameter_vector_.size(); i++) {
177 RCP<Epetra_Map> local_map = rcp(new Epetra_LocalMap(static_cast<int>(parameter_vector_[i].size()), 0, map_x_->Comm())) ;
178 p_map_.push_back(local_map);
179
180 RCP<Epetra_Vector> ep_vec = rcp(new Epetra_Vector(*local_map));
181 ep_vec->PutScalar(0.0);
182
183 for (unsigned int j=0; j < parameter_vector_[i].size(); j++)
184 (*ep_vec)[j] = parameter_vector_[i][j].baseValue;
185
186 p_init_.push_back(ep_vec);
187 }
188
189 // Initialize the epetra operator factory
190 epetraOperatorFactory_ = Teuchos::rcp(new EpetraLOF_EOpFactory(lof));
191}
192
194initializeParameterVector(const std::vector<Teuchos::RCP<Teuchos::Array<std::string> > >& p_names,
195 const std::vector<Teuchos::RCP<Teuchos::Array<double> > >& p_values,
196 const Teuchos::RCP<panzer::ParamLib>& parameter_library)
197{
198 parameter_vector_.resize(p_names.size());
199 is_distributed_parameter_.resize(p_names.size(),false);
200 for (std::vector<Teuchos::RCP<Teuchos::Array<std::string> > >::size_type p = 0;
201 p < p_names.size(); ++p) {
202
203 // register all the scalar parameters, setting initial
204 for(int i=0;i<p_names[p]->size();i++)
205 registerScalarParameter((*p_names[p])[i],*parameter_library,(*p_values[p])[i]);
206
207 parameter_library->fillVector<panzer::Traits::Residual>(*(p_names[p]), parameter_vector_[p]);
208
209 for(int i=0;i<p_names[p]->size();i++) {
210 parameter_vector_[p][i].baseValue = (*p_values[p])[i];
211 parameter_vector_[p][i].family->setRealValueForAllTypes((*p_values[p])[i]);
212 }
213 }
214}
215
216// Post-Construction methods to add parameters and responses
217
219addDistributedParameter(const std::string name,
220 const Teuchos::RCP<Epetra_Map>& global_map,
221 const Teuchos::RCP<Epetra_Import>& importer,
222 const Teuchos::RCP<Epetra_Vector>& ghosted_vector)
223{
224 // Will push_back a new parameter entry
225 int index = static_cast<int>(p_map_.size());
226
227 p_map_.push_back(global_map);
228 Teuchos::RCP<Epetra_Vector> ep_vec = Teuchos::rcp(new Epetra_Vector(*global_map));
229 ep_vec->PutScalar(0.0);
230 p_init_.push_back(ep_vec);
231
232 Teuchos::RCP<Teuchos::Array<std::string> > tmp_names =
233 Teuchos::rcp(new Teuchos::Array<std::string>);
234 tmp_names->push_back(name);
235 p_names_.push_back(tmp_names);
236
237 is_distributed_parameter_.push_back(true);
238
239 // NOTE: we do not add this parameter to the sacado parameter
240 // library in the global data object. That library is for scalars.
241 // We will need to do something different if we need sensitivities
242 // wrt distributed parameters.
243
244 distributed_parameter_container_.push_back(std::make_tuple(name,index,importer,ghosted_vector));
245
246 return index;
247}
248
249// Overridden from EpetraExt::ModelEvaluator
250
251Teuchos::RCP<const Epetra_Map>
253{
254 return map_x_;
255}
256
257Teuchos::RCP<const Epetra_Map>
259{
260 return map_x_;
261}
262
263Teuchos::RCP<const Epetra_Vector>
265{
266 return x0_;
267}
268
269Teuchos::RCP<const Epetra_Vector>
271{
272 return x_dot_init_;
273}
274
275double
277{
278 return t_init_;
279}
280
281Teuchos::RCP<Epetra_Operator>
283{
284 return epetraOperatorFactory_->create();
285}
286
287Teuchos::RCP<const Epetra_Map>
289{
290 return p_map_[l];
291}
292
293Teuchos::RCP<const Teuchos::Array<std::string> >
295{
296 return p_names_[l];
297}
298
299Teuchos::RCP<const Epetra_Vector>
301{
302 return p_init_[l];
303}
304
305Teuchos::RCP<const Epetra_Map>
307{
308 return g_map_[l];
309}
310
313{
314 InArgsSetup inArgs;
315 inArgs.setModelEvalDescription(this->description());
316 inArgs.setSupports(IN_ARG_x,true);
317 if (build_transient_support_) {
318 inArgs.setSupports(IN_ARG_x_dot,true);
319 inArgs.setSupports(IN_ARG_t,true);
320 inArgs.setSupports(IN_ARG_alpha,true);
321 inArgs.setSupports(IN_ARG_beta,true);
322 }
323 inArgs.set_Np(p_init_.size());
324
325 return inArgs;
326}
327
330{
331 OutArgsSetup outArgs;
332 outArgs.setModelEvalDescription(this->description());
333 outArgs.set_Np_Ng(p_init_.size(), g_map_.size());
334 outArgs.setSupports(OUT_ARG_f,true);
335 outArgs.setSupports(OUT_ARG_W,true);
336 outArgs.set_W_properties(
337 DerivativeProperties(
338 DERIV_LINEARITY_NONCONST
339 ,DERIV_RANK_FULL
340 ,true // supportsAdjoint
341 )
342 );
343
344 // add in df/dp (if appropriate)
345 for(std::size_t i=0;i<p_init_.size();i++) {
346 if(!is_distributed_parameter_[i])
347 outArgs.setSupports(OUT_ARG_DfDp,i,EpetraExt::ModelEvaluator::DerivativeSupport(DERIV_MV_BY_COL));
348 }
349
350 // add in dg/dx (if appropriate)
351 for(std::size_t i=0;i<g_names_.size();i++) {
352 typedef panzer::Traits::Jacobian RespEvalT;
353
354 // check dg/dx and add it in if appropriate
355 Teuchos::RCP<panzer::ResponseBase> respJacBase = responseLibrary_->getResponse<RespEvalT>(g_names_[i]);
356 if(respJacBase!=Teuchos::null) {
357 // cast is guranteed to succeed because of check in addResponse
358 Teuchos::RCP<panzer::ResponseMESupportBase<RespEvalT> > resp = Teuchos::rcp_dynamic_cast<panzer::ResponseMESupportBase<RespEvalT> >(respJacBase);
359
360 // class must supoprt a derivative
361 if(resp->supportsDerivative())
362 outArgs.setSupports(OUT_ARG_DgDx,i,DerivativeSupport(DERIV_TRANS_MV_BY_ROW));
363 //outArgs.setSupports(OUT_ARG_DgDx,i,DerivativeSupport(DERIV_LINEAR_OP));
364 }
365 }
366
367 return outArgs;
368}
369
371applyDirichletBCs(const Teuchos::RCP<Thyra::VectorBase<double> > & x,
372 const Teuchos::RCP<Thyra::VectorBase<double> > & f) const
373{
374 using Teuchos::RCP;
375 using Teuchos::ArrayRCP;
376 using Teuchos::Array;
377 //using Teuchos::tuple;
378 using Teuchos::rcp_dynamic_cast;
379
380 // if neccessary build a ghosted container
381 if(Teuchos::is_null(ghostedContainer_)) {
382 ghostedContainer_ = lof_->buildGhostedLinearObjContainer();
383 lof_->initializeGhostedContainer(panzer::LinearObjContainer::X |
384 panzer::LinearObjContainer::F,*ghostedContainer_);
385 }
386
388 ae_inargs.container_ = lof_->buildLinearObjContainer(); // we use a new global container
389 ae_inargs.ghostedContainer_ = ghostedContainer_; // we can reuse the ghosted container
390 ae_inargs.alpha = 0.0;
391 ae_inargs.beta = 1.0;
392 ae_inargs.evaluate_transient_terms = false;
393
394 // this is the tempory target
395 lof_->initializeContainer(panzer::LinearObjContainer::F,*ae_inargs.container_);
396
397 // here we are building a container, this operation is fast, simply allocating a struct
398 const RCP<panzer::ThyraObjContainer<double> > thGlobalContainer =
399 Teuchos::rcp_dynamic_cast<panzer::ThyraObjContainer<double> >(ae_inargs.container_,true);
400
401 TEUCHOS_ASSERT(!Teuchos::is_null(thGlobalContainer));
402
403 // Ghosted container objects are zeroed out below only if needed for
404 // a particular calculation. This makes it more efficient than
405 // zeroing out all objects in the container here.
406 const RCP<panzer::ThyraObjContainer<double> > thGhostedContainer =
407 Teuchos::rcp_dynamic_cast<panzer::ThyraObjContainer<double> >(ae_inargs.ghostedContainer_,true);
408 TEUCHOS_ASSERT(!Teuchos::is_null(thGhostedContainer));
409 Thyra::assign(thGhostedContainer->get_f_th().ptr(),0.0);
410
411 // Set the solution vector (currently all targets require solution).
412 // In the future we may move these into the individual cases below.
413 // A very subtle (and fragile) point: A non-null pointer in global
414 // container triggers export operations during fill. Also, the
415 // introduction of the container is forcing us to cast away const on
416 // arguments that should be const. Another reason to redesign
417 // LinearObjContainer layers.
418 thGlobalContainer->set_x_th(x);
419
420 // evaluate dirichlet boundary conditions
421 RCP<panzer::LinearObjContainer> counter = ae_tm_.getAsObject<panzer::Traits::Residual>()->evaluateOnlyDirichletBCs(ae_inargs);
422
423 // allocate the result container
424 RCP<panzer::LinearObjContainer> result = lof_->buildLinearObjContainer(); // we use a new global container
425
426 // stuff the evaluate boundary conditions into the f spot of the counter ... the x is already filled
427 Teuchos::rcp_dynamic_cast<panzer::ThyraObjContainer<double> >(counter)->set_f_th(thGlobalContainer->get_f_th());
428
429 // stuff the vector that needs applied dirichlet conditions in the the f spot of the result LOC
430 Teuchos::rcp_dynamic_cast<panzer::ThyraObjContainer<double> >(result)->set_f_th(f);
431
432 // use the linear object factory to apply the result
433 lof_->applyDirichletBCs(*counter,*result);
434}
435
437 const OutArgs& outArgs ) const
438{
439 using Teuchos::RCP;
440 using Teuchos::rcp_dynamic_cast;
441
442 evalModel_basic(inArgs,outArgs);
443}
444
446 const OutArgs& outArgs ) const
447{
448 using Teuchos::RCP;
449 using Teuchos::rcp_dynamic_cast;
450
451 // Transient or steady-state evaluation is determined by the x_dot
452 // vector. If this RCP is null, then we are doing a steady-state
453 // fill.
454 bool has_x_dot = false;
455 if (inArgs.supports(EpetraExt::ModelEvaluator::IN_ARG_x_dot ))
456 has_x_dot = nonnull(inArgs.get_x_dot());
457
458 // Make sure construction built in transient support
459 TEUCHOS_TEST_FOR_EXCEPTION(has_x_dot && !build_transient_support_, std::runtime_error,
460 "ModelEvaluator was not built with transient support enabled!");
461
462 //
463 // Get the output arguments
464 //
465 const RCP<Epetra_Vector> f_out = outArgs.get_f();
466 const RCP<Epetra_Operator> W_out = outArgs.get_W();
467 bool requiredResponses = required_basic_g(outArgs);
468 bool requiredSensitivities = required_basic_dfdp(outArgs);
469
470 // see if the user wants us to do anything
471 if(Teuchos::is_null(f_out) && Teuchos::is_null(W_out) &&
472 !requiredResponses && !requiredSensitivities) {
473 return;
474 }
475
476 // the user requested work from this method
477 // keep on moving
478
479 // if neccessary build a ghosted container
480 if(Teuchos::is_null(ghostedContainer_)) {
481 ghostedContainer_ = lof_->buildGhostedLinearObjContainer();
482 lof_->initializeGhostedContainer(panzer::LinearObjContainer::X |
485 panzer::LinearObjContainer::Mat, *ghostedContainer_);
486 }
487
488 //
489 // Get the input arguments
490 //
491 const RCP<const Epetra_Vector> x = inArgs.get_x();
492 RCP<const Epetra_Vector> x_dot;
493
495 ae_inargs.container_ = lof_->buildLinearObjContainer(); // we use a new global container
496 ae_inargs.ghostedContainer_ = ghostedContainer_; // we can reuse the ghosted container
497 ae_inargs.alpha = 0.0;
498 ae_inargs.beta = 1.0;
499 ae_inargs.evaluate_transient_terms = false;
500 if (build_transient_support_) {
501 x_dot = inArgs.get_x_dot();
502 ae_inargs.alpha = inArgs.get_alpha();
503 ae_inargs.beta = inArgs.get_beta();
504 ae_inargs.time = inArgs.get_t();
505 ae_inargs.evaluate_transient_terms = true;
506 }
507
508 // handle application of the one time dirichlet beta in the
509 // assembly engine. Note that this has to be set explicitly
510 // each time because this badly breaks encapsulation. Essentially
511 // we must work around the model evaluator abstraction!
512 ae_inargs.apply_dirichlet_beta = false;
513 if(oneTimeDirichletBeta_on_) {
514 ae_inargs.dirichlet_beta = oneTimeDirichletBeta_;
515 ae_inargs.apply_dirichlet_beta = true;
516
517 oneTimeDirichletBeta_on_ = false;
518 }
519
520 // Set locally replicated scalar input parameters
521 for (int i=0; i<inArgs.Np(); i++) {
522 Teuchos::RCP<const Epetra_Vector> p = inArgs.get_p(i);
523 if ( nonnull(p) && !is_distributed_parameter_[i]) {
524 for (unsigned int j=0; j < parameter_vector_[i].size(); j++) {
525 parameter_vector_[i][j].baseValue = (*p)[j];
526 }
527 }
528 }
529
530 for (Teuchos::Array<panzer::ParamVec>::size_type i=0; i < parameter_vector_.size(); i++) {
531 for (unsigned int j=0; j < parameter_vector_[i].size(); j++) {
532 parameter_vector_[i][j].family->setRealValueForAllTypes(parameter_vector_[i][j].baseValue);
533 }
534 }
535
536 // Perform global to ghost and set distributed parameters
537 for (std::vector<std::tuple<std::string,int,Teuchos::RCP<Epetra_Import>,Teuchos::RCP<Epetra_Vector> > >::const_iterator i =
538 distributed_parameter_container_.begin(); i != distributed_parameter_container_.end(); ++i) {
539 // do export if parameter exists in inArgs
540 Teuchos::RCP<const Epetra_Vector> global_vec = inArgs.get_p(std::get<1>(*i));
541 if (nonnull(global_vec)) {
542 // Only import if the importer is nonnull
543 Teuchos::RCP<Epetra_Import> importer = std::get<2>(*i);
544 if (nonnull(importer))
545 std::get<3>(*i)->Import(*global_vec,*importer,Insert);
546 }
547
548 // set in ae_inargs_ string lookup container
549 Teuchos::RCP<panzer::EpetraLinearObjContainer> container =
550 Teuchos::rcp(new panzer::EpetraLinearObjContainer(p_map_[std::get<1>(*i)],p_map_[std::get<1>(*i)]));
551 container->set_x(std::get<3>(*i));
552 ae_inargs.addGlobalEvaluationData(std::get<0>(*i),container);
553 }
554
555 // here we are building a container, this operation is fast, simply allocating a struct
556 const RCP<panzer::EpetraLinearObjContainer> epGlobalContainer =
557 Teuchos::rcp_dynamic_cast<panzer::EpetraLinearObjContainer>(ae_inargs.container_);
558
559 TEUCHOS_ASSERT(!Teuchos::is_null(epGlobalContainer));
560
561 // Ghosted container objects are zeroed out below only if needed for
562 // a particular calculation. This makes it more efficient thatn
563 // zeroing out all objects in the container here.
564 const RCP<panzer::EpetraLinearObjContainer> epGhostedContainer =
565 Teuchos::rcp_dynamic_cast<panzer::EpetraLinearObjContainer>(ae_inargs.ghostedContainer_);
566
567 // Set the solution vector (currently all targets require solution).
568 // In the future we may move these into the individual cases below.
569 // A very subtle (and fragile) point: A non-null pointer in global
570 // container triggers export operations during fill. Also, the
571 // introduction of the container is forcing us to cast away const on
572 // arguments that should be const. Another reason to redesign
573 // LinearObjContainer layers.
574 epGlobalContainer->set_x(Teuchos::rcp_const_cast<Epetra_Vector>(x));
575 if (has_x_dot)
576 epGlobalContainer->set_dxdt(Teuchos::rcp_const_cast<Epetra_Vector>(x_dot));
577
578 if (!Teuchos::is_null(f_out) && !Teuchos::is_null(W_out)) {
579
580 PANZER_FUNC_TIME_MONITOR("panzer::ModelEvaluator::evalModel(f and J)");
581
582 // Set the targets
583 epGlobalContainer->set_f(f_out);
584 epGlobalContainer->set_A(Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(W_out));
585
586 // Zero values in ghosted container objects
587 epGhostedContainer->get_f()->PutScalar(0.0);
588 epGhostedContainer->get_A()->PutScalar(0.0);
589
590 ae_tm_.getAsObject<panzer::Traits::Jacobian>()->evaluate(ae_inargs);
591 }
592 else if(!Teuchos::is_null(f_out) && Teuchos::is_null(W_out)) {
593
594 PANZER_FUNC_TIME_MONITOR("panzer::ModelEvaluator::evalModel(f)");
595
596 epGlobalContainer->set_f(f_out);
597
598 // Zero values in ghosted container objects
599 epGhostedContainer->get_f()->PutScalar(0.0);
600
601 ae_tm_.getAsObject<panzer::Traits::Residual>()->evaluate(ae_inargs);
602
603 f_out->Update(1.0, *(epGlobalContainer->get_f()), 0.0);
604 }
605 else if(Teuchos::is_null(f_out) && !Teuchos::is_null(W_out)) {
606
607 PANZER_FUNC_TIME_MONITOR("panzer::ModelEvaluator::evalModel(J)");
608
609 // this dummy nonsense is needed only for scattering dirichlet conditions
610 if (Teuchos::is_null(dummy_f_))
611 dummy_f_ = Teuchos::rcp(new Epetra_Vector(*(this->get_f_map())));
612 epGlobalContainer->set_f(dummy_f_);
613 epGlobalContainer->set_A(Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(W_out));
614
615 // Zero values in ghosted container objects
616 epGhostedContainer->get_A()->PutScalar(0.0);
617
618 ae_tm_.getAsObject<panzer::Traits::Jacobian>()->evaluate(ae_inargs);
619 }
620 // HACK: set A to null before calling responses to avoid touching the
621 // the Jacobian after it has been properly assembled. Should be fixed
622 // by using a modified version of ae_inargs instead.
623 epGlobalContainer->set_A(Teuchos::null);
624
625 // evaluate responses...uses the stored assembly arguments and containers
626 if(requiredResponses) {
627 evalModel_basic_g(ae_inargs,inArgs,outArgs);
628
629 // evaluate response derivatives
630 if(required_basic_dgdx(outArgs))
631 evalModel_basic_dgdx(ae_inargs,inArgs,outArgs);
632 }
633
634 if(required_basic_dfdp(outArgs))
635 evalModel_basic_dfdp(ae_inargs,inArgs,outArgs);
636
637 // Holding a rcp to f produces a seg fault in Rythmos when the next
638 // f comes in and the resulting dtor is called. Need to discuss
639 // with Ross. Clearing all references here works!
640
641 epGlobalContainer->set_x(Teuchos::null);
642 epGlobalContainer->set_dxdt(Teuchos::null);
643 epGlobalContainer->set_f(Teuchos::null);
644 epGlobalContainer->set_A(Teuchos::null);
645
646 // forget previous containers
647 ae_inargs.container_ = Teuchos::null;
648 ae_inargs.ghostedContainer_ = Teuchos::null;
649}
650
651void
653evalModel_basic_g(AssemblyEngineInArgs ae_inargs, const InArgs& /* inArgs */, const OutArgs& outArgs) const
654{
655 // optional sanity check
656 // TEUCHOS_ASSERT(required_basic_g(outArgs));
657
658 for(std::size_t i=0;i<g_names_.size();i++) {
659 Teuchos::RCP<Epetra_Vector> vec = outArgs.get_g(i);
660 if(vec!=Teuchos::null) {
661 std::string responseName = g_names_[i];
662 Teuchos::RCP<panzer::ResponseMESupportBase<panzer::Traits::Residual> > resp
663 = Teuchos::rcp_dynamic_cast<panzer::ResponseMESupportBase<panzer::Traits::Residual> >(responseLibrary_->getResponse<panzer::Traits::Residual>(responseName));
664 resp->setVector(vec);
665 }
666 }
667
668 // evaluator responses
669 responseLibrary_->addResponsesToInArgs<panzer::Traits::Residual>(ae_inargs);
670 responseLibrary_->evaluate<panzer::Traits::Residual>(ae_inargs);
671}
672
673void
675evalModel_basic_dgdx(AssemblyEngineInArgs ae_inargs, const InArgs& /* inArgs */, const OutArgs& outArgs) const
676{
677 // optional sanity check
678 TEUCHOS_ASSERT(required_basic_dgdx(outArgs));
679
680 for(std::size_t i=0;i<g_names_.size();i++) {
681 // get "Vector" out of derivative, if its something else, throw an exception
682 EpetraExt::ModelEvaluator::Derivative deriv = outArgs.get_DgDx(i);
683 if(deriv.isEmpty())
684 continue;
685
686 Teuchos::RCP<Epetra_MultiVector> vec = deriv.getMultiVector();
687
688 if(vec!=Teuchos::null) {
689 std::string responseName = g_names_[i];
690 Teuchos::RCP<panzer::ResponseMESupportBase<panzer::Traits::Jacobian> > resp
691 = Teuchos::rcp_dynamic_cast<panzer::ResponseMESupportBase<panzer::Traits::Jacobian> >(responseLibrary_->getResponse<panzer::Traits::Jacobian>(responseName));
692 resp->setDerivative(vec);
693 }
694 }
695
696 // evaluator responses
697 responseLibrary_->addResponsesToInArgs<panzer::Traits::Jacobian>(ae_inargs);
698 responseLibrary_->evaluate<panzer::Traits::Jacobian>(ae_inargs);
699}
700
701void
703evalModel_basic_dfdp(AssemblyEngineInArgs ae_inargs, const InArgs& /* inArgs */, const OutArgs& outArgs) const
704{
705 using Teuchos::RCP;
706 using Teuchos::rcp_dynamic_cast;
707
708 TEUCHOS_ASSERT(required_basic_dfdp(outArgs));
709
710 // dynamic cast to blocked LOF for now
711 RCP<const Thyra::VectorSpaceBase<double> > glblVS = rcp_dynamic_cast<const ThyraObjFactory<double> >(lof_,true)->getThyraRangeSpace();;
712
713 std::vector<std::string> activeParameters;
714
715 // fill parameter vector containers
716 int totalParameterCount = 0;
717 for(Teuchos::Array<panzer::ParamVec>::size_type i=0; i < parameter_vector_.size(); i++) {
718 // have derivatives been requested?
719 EpetraExt::ModelEvaluator::Derivative deriv = outArgs.get_DfDp(i);
720 if(deriv.isEmpty())
721 continue;
722
723 // grab multivector, make sure its the right dimension
724 Teuchos::RCP<Epetra_MultiVector> mVec = deriv.getMultiVector();
725 TEUCHOS_ASSERT(mVec->NumVectors()==int(parameter_vector_[i].size()));
726
727 for (unsigned int j=0; j < parameter_vector_[i].size(); j++) {
728
729 // build containers for each vector
730 RCP<LOCPair_GlobalEvaluationData> loc_pair = Teuchos::rcp(new LOCPair_GlobalEvaluationData(lof_,LinearObjContainer::F));
731 RCP<LinearObjContainer> globalContainer = loc_pair->getGlobalLOC();
732
733 // stuff target vector into global container
734 RCP<Epetra_Vector> vec = Teuchos::rcpFromRef(*(*mVec)(j));
735 RCP<panzer::ThyraObjContainer<double> > thGlobalContainer =
736 Teuchos::rcp_dynamic_cast<panzer::ThyraObjContainer<double> >(globalContainer);
737 thGlobalContainer->set_f_th(Thyra::create_Vector(vec,glblVS));
738
739 // add container into in args object
740 std::string name = "PARAMETER_SENSITIVIES: "+(*p_names_[i])[j];
741 ae_inargs.addGlobalEvaluationData(name,loc_pair->getGhostedLOC());
742 ae_inargs.addGlobalEvaluationData(name+"_pair",loc_pair);
743
744 activeParameters.push_back(name);
745 totalParameterCount++;
746 }
747 }
748
749 // this communicates to the scatter evaluators so that the appropriate parameters are scattered
750 RCP<GlobalEvaluationData> ged_activeParameters = Teuchos::rcp(new ParameterList_GlobalEvaluationData(activeParameters));
751 ae_inargs.addGlobalEvaluationData("PARAMETER_NAMES",ged_activeParameters);
752
753 int paramIndex = 0;
754 for (Teuchos::Array<panzer::ParamVec>::size_type i=0; i < parameter_vector_.size(); i++) {
755 // don't modify the parameter if its not needed
756 EpetraExt::ModelEvaluator::Derivative deriv = outArgs.get_DfDp(i);
757 if(deriv.isEmpty()) {
758 // reinitialize values that should not have sensitivities computed (this is a precaution)
759 for (unsigned int j=0; j < parameter_vector_[i].size(); j++) {
760 Traits::FadType p = Traits::FadType(totalParameterCount, parameter_vector_[i][j].baseValue);
761 parameter_vector_[i][j].family->setValue<panzer::Traits::Tangent>(p);
762 }
763 continue;
764 }
765 else {
766 // loop over each parameter in the vector, initializing the AD type
767 for (unsigned int j=0; j < parameter_vector_[i].size(); j++) {
768 Traits::FadType p = Traits::FadType(totalParameterCount, parameter_vector_[i][j].baseValue);
769 p.fastAccessDx(paramIndex) = 1.0;
770 parameter_vector_[i][j].family->setValue<panzer::Traits::Tangent>(p);
771 paramIndex++;
772 }
773 }
774 }
775
776 // make sure that the total parameter count and the total parameter index match up
777 TEUCHOS_ASSERT(paramIndex==totalParameterCount);
778
779 if(totalParameterCount>0) {
780 ae_tm_.getAsObject<panzer::Traits::Tangent>()->evaluate(ae_inargs);
781 }
782}
783
784bool panzer::ModelEvaluator_Epetra::required_basic_g(const OutArgs & outArgs) const
785{
786 // determine if any of the outArgs are not null!
787 bool activeGArgs = false;
788 for(int i=0;i<outArgs.Ng();i++)
789 activeGArgs |= (outArgs.get_g(i)!=Teuchos::null);
790
791 return activeGArgs | required_basic_dgdx(outArgs);
792}
793
794bool panzer::ModelEvaluator_Epetra::required_basic_dgdx(const OutArgs & outArgs) const
795{
796 // determine if any of the outArgs are not null!
797 bool activeGArgs = false;
798 for(int i=0;i<outArgs.Ng();i++) {
799 // no derivatives are supported
800 if(outArgs.supports(OUT_ARG_DgDx,i).none())
801 continue;
802
803 // this is basically a redundant computation
804 activeGArgs |= (!outArgs.get_DgDx(i).isEmpty());
805 }
806
807 return activeGArgs;
808}
809
810bool panzer::ModelEvaluator_Epetra::required_basic_dfdp(const OutArgs & outArgs) const
811{
812 // determine if any of the outArgs are not null!
813 bool activeFPArgs = false;
814 for(int i=0;i<outArgs.Np();i++) {
815 // no derivatives are supported
816 if(outArgs.supports(OUT_ARG_DfDp,i).none())
817 continue;
818
819 // this is basically a redundant computation
820 activeFPArgs |= (!outArgs.get_DfDp(i).isEmpty());
821 }
822
823 return activeFPArgs;
824}
825
827{
828 t_init_ = t;
829}
830
833{
834
835 using Teuchos::rcpFromRef;
836 using Teuchos::rcp_dynamic_cast;
837 using Teuchos::RCP;
838 using Teuchos::ArrayView;
839 using Teuchos::rcp_dynamic_cast;
840
841 const int numVecs = x.NumVectors();
842
843 TEUCHOS_TEST_FOR_EXCEPTION(numVecs != 1, std::runtime_error,
844 "epetraToThyra does not work with MV dimension != 1");
845
846 const RCP<const Thyra::ProductVectorBase<double> > prodThyraVec =
847 Thyra::castOrCreateProductVectorBase(rcpFromRef(thyraVec));
848
849 const Teuchos::ArrayView<double> epetraData(x[0], x.Map().NumMyElements());
850 // NOTE: See above!
851
852 std::size_t offset = 0;
853 const int numBlocks = prodThyraVec->productSpace()->numBlocks();
854 for (int b = 0; b < numBlocks; ++b) {
855 const RCP<const Thyra::VectorBase<double> > vec_b = prodThyraVec->getVectorBlock(b);
856 const RCP<const Thyra::SpmdVectorBase<double> > spmd_b =
857 rcp_dynamic_cast<const Thyra::SpmdVectorBase<double> >(vec_b, true);
858
859 Teuchos::ArrayRCP<const double> thyraData;
860 spmd_b->getLocalData(Teuchos::ptrFromRef(thyraData));
861
862 for (Teuchos::ArrayRCP<const double>::size_type i=0; i < thyraData.size(); ++i) {
863 epetraData[i+offset] = thyraData[i];
864 }
865 offset += thyraData.size();
866 }
867
868}
869
870
872copyEpetraIntoThyra(const Epetra_MultiVector& x, const Teuchos::Ptr<Thyra::VectorBase<double> > &thyraVec) const
873{
874 using Teuchos::RCP;
875 using Teuchos::ArrayView;
876 using Teuchos::rcpFromPtr;
877 using Teuchos::rcp_dynamic_cast;
878
879 const int numVecs = x.NumVectors();
880
881 TEUCHOS_TEST_FOR_EXCEPTION(numVecs != 1, std::runtime_error,
882 "ModelEvaluator_Epetra::copyEpetraToThyra does not work with MV dimension != 1");
883
884 const RCP<Thyra::ProductVectorBase<double> > prodThyraVec =
885 Thyra::castOrCreateNonconstProductVectorBase(rcpFromPtr(thyraVec));
886
887 const Teuchos::ArrayView<const double> epetraData(x[0], x.Map().NumMyElements());
888 // NOTE: I tried using Epetra_MultiVector::operator()(int) to return an
889 // Epetra_Vector object but it has a defect when Reset(...) is called which
890 // results in a memory access error (see bug 4700).
891
892 std::size_t offset = 0;
893 const int numBlocks = prodThyraVec->productSpace()->numBlocks();
894 for (int b = 0; b < numBlocks; ++b) {
895 const RCP<Thyra::VectorBase<double> > vec_b = prodThyraVec->getNonconstVectorBlock(b);
896 const RCP<Thyra::SpmdVectorBase<double> > spmd_b =
897 rcp_dynamic_cast<Thyra::SpmdVectorBase<double> >(vec_b, true);
898
899 Teuchos::ArrayRCP<double> thyraData;
900 spmd_b->getNonconstLocalData(Teuchos::ptrFromRef(thyraData));
901
902 for (Teuchos::ArrayRCP<double>::size_type i=0; i < thyraData.size(); ++i) {
903 thyraData[i] = epetraData[i+offset];
904 }
905 offset += thyraData.size();
906 }
907
908}
909
910
911Teuchos::RCP<panzer::ModelEvaluator_Epetra>
912panzer::buildEpetraME(const Teuchos::RCP<panzer::FieldManagerBuilder>& fmb,
913 const Teuchos::RCP<panzer::ResponseLibrary<panzer::Traits> >& rLibrary,
914 const Teuchos::RCP<panzer::LinearObjFactory<panzer::Traits> >& lof,
915 const std::vector<Teuchos::RCP<Teuchos::Array<std::string> > >& p_names,
916 const std::vector<Teuchos::RCP<Teuchos::Array<double> > >& p_values,
917 const Teuchos::RCP<panzer::GlobalData>& global_data,
918 bool build_transient_support)
919{
920 using Teuchos::RCP;
921 using Teuchos::rcp_dynamic_cast;
922 using Teuchos::rcp;
923
924 std::stringstream ss;
925 ss << "panzer::buildEpetraME: Linear object factory is incorrect type. Should be one of: ";
926
927 RCP<BlockedEpetraLinearObjFactory<panzer::Traits,int> > ep_lof
928 = rcp_dynamic_cast<BlockedEpetraLinearObjFactory<panzer::Traits,int> >(lof);
929
930 // if you can, build from an epetra linear object factory
931 if(ep_lof!=Teuchos::null)
932 return rcp(new ModelEvaluator_Epetra(fmb,rLibrary,ep_lof,p_names,p_values,global_data,build_transient_support));
933
934 ss << "\"BlockedEpetraLinearObjFactory\", ";
935
936 TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,ss.str());
937}
938
940setOneTimeDirichletBeta(const double & beta) const
941{
942 oneTimeDirichletBeta_on_ = true;
943 oneTimeDirichletBeta_ = beta;
944}
Insert
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.
Teuchos::RCP< Epetra_MultiVector > getMultiVector() const
int NumMyElements() const
const Epetra_BlockMap & Map() const
int NumVectors() const
void addGlobalEvaluationData(const std::string &key, const Teuchos::RCP< GlobalEvaluationData > &ged)
Teuchos::RCP< panzer::LinearObjContainer > ghostedContainer_
Teuchos::RCP< panzer::LinearObjContainer > container_
virtual const Teuchos::RCP< Epetra_Map > getMap(int i) const
get the map from the matrix
Teuchos::RCP< const Epetra_Map > get_f_map() const
Teuchos::RCP< const Epetra_Map > get_x_map() const
Teuchos::RCP< const Epetra_Vector > get_x_init() const
Teuchos::RCP< const Epetra_Vector > get_x_dot_init() const
Teuchos::RCP< Epetra_Operator > create_W() const
Teuchos::RCP< const Epetra_Map > get_g_map(int l) const
void evalModel(const InArgs &inArgs, const OutArgs &outArgs) const
Teuchos::RCP< const Teuchos::Array< std::string > > get_p_names(int l) const
void set_t_init(double t)
Set initial time value.
ModelEvaluator_Epetra(const Teuchos::RCP< panzer::FieldManagerBuilder > &fmb, const Teuchos::RCP< panzer::ResponseLibrary< panzer::Traits > > &rLibrary, const Teuchos::RCP< panzer::LinearObjFactory< panzer::Traits > > &lof, const std::vector< Teuchos::RCP< Teuchos::Array< std::string > > > &p_names, const std::vector< Teuchos::RCP< Teuchos::Array< double > > > &p_values, const Teuchos::RCP< panzer::GlobalData > &global_data, bool build_transient_support)
Teuchos::RCP< const Epetra_Map > get_p_map(int l) const
Teuchos::RCP< const Epetra_Vector > get_p_init(int l) const
void evalModel_basic_dgdx(AssemblyEngineInArgs ae_inargs, const InArgs &inArgs, const OutArgs &outArgs) const
void initializeEpetraObjs(panzer::BlockedEpetraLinearObjFactory< panzer::Traits, int > &lof)
void setOneTimeDirichletBeta(const double &beta) const
bool required_basic_dgdx(const OutArgs &outArgs) const
Are their required responses in the out args? DgDx.
void evalModel_basic_dfdp(AssemblyEngineInArgs ae_inargs, const InArgs &inArgs, const OutArgs &outArgs) const
bool required_basic_g(const OutArgs &outArgs) const
Are their required responses in the out args? g and DgDx.
int addDistributedParameter(const std::string name, const Teuchos::RCP< Epetra_Map > &global_map, const Teuchos::RCP< Epetra_Import > &importer, const Teuchos::RCP< Epetra_Vector > &ghosted_vector)
void evalModel_basic_g(AssemblyEngineInArgs ae_inargs, const InArgs &inArgs, const OutArgs &outArgs) const
bool required_basic_dfdp(const OutArgs &outArgs) const
Are derivatives of the residual with respect to the parameters in the out args? DfDp.
void evalModel_basic(const InArgs &inArgs, const OutArgs &outArgs) const
for evaluation and handling of normal quantities, x,f,W, etc
void initializeParameterVector(const std::vector< Teuchos::RCP< Teuchos::Array< std::string > > > &p_names, const std::vector< Teuchos::RCP< Teuchos::Array< double > > > &p_values, const Teuchos::RCP< panzer::ParamLib > &parameter_library)
void copyEpetraIntoThyra(const Epetra_MultiVector &x, const Teuchos::Ptr< Thyra::VectorBase< double > > &thyraVec) const
panzer::AssemblyEngine_TemplateManager< panzer::Traits > ae_tm_
void applyDirichletBCs(const Teuchos::RCP< Thyra::VectorBase< double > > &x, const Teuchos::RCP< Thyra::VectorBase< double > > &f) const
std::vector< Teuchos::RCP< Teuchos::Array< std::string > > > p_names_
void copyThyraIntoEpetra(const Thyra::VectorBase< double > &thyraVec, Epetra_MultiVector &x) const
Teuchos::RCP< ModelEvaluator_Epetra > buildEpetraME(const Teuchos::RCP< FieldManagerBuilder > &fmb, const Teuchos::RCP< ResponseLibrary< panzer::Traits > > &rLibrary, const Teuchos::RCP< LinearObjFactory< panzer::Traits > > &lof, const std::vector< Teuchos::RCP< Teuchos::Array< std::string > > > &p_names, const std::vector< Teuchos::RCP< Teuchos::Array< double > > > &p_values, const Teuchos::RCP< panzer::GlobalData > &global_data, bool build_transient_support)
void registerScalarParameter(const std::string name, panzer::ParamLib &pl, double realValue)
PANZER_FADTYPE FadType