Panzer  Version of the Day
Panzer_STK_ModelEvaluatorFactory_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_STK_MODEL_EVALUATOR_FACTORY_T_HPP
44 #define PANZER_STK_MODEL_EVALUATOR_FACTORY_T_HPP
45 
46 #include "Thyra_ModelEvaluator.hpp"
47 #include "Teuchos_Assert.hpp"
48 #include "Teuchos_as.hpp"
49 #include "Teuchos_DefaultMpiComm.hpp"
50 #include "Teuchos_AbstractFactoryStd.hpp"
51 
52 #include "PanzerAdaptersSTK_config.hpp"
53 #include "Panzer_Traits.hpp"
54 #include "Panzer_GlobalData.hpp"
55 #include "Panzer_BC.hpp"
57 #include "Panzer_BasisIRLayout.hpp"
58 #include "Panzer_DOFManager.hpp"
63 #include "Panzer_TpetraLinearObjFactory.hpp"
77 
78 #include "Panzer_STK_Interface.hpp"
88 #include "Panzer_STK_Utilities.hpp"
98 
99 #include <vector>
100 #include <iostream>
101 #include <fstream>
102 #include <cstdlib> // for std::getenv
103 
104 // Piro solver objects
105 #include "Thyra_EpetraModelEvaluator.hpp"
106 #include "Piro_ConfigDefs.hpp"
107 #include "Piro_NOXSolver.hpp"
108 #include "Piro_LOCASolver.hpp"
109 #include "Piro_RythmosSolver.hpp"
110 
111 #include <Panzer_NodeType.hpp>
112 
113 namespace panzer_stk {
114 
115  template<typename ScalarT>
116  void ModelEvaluatorFactory<ScalarT>::setParameterList(Teuchos::RCP<Teuchos::ParameterList> const& paramList)
117  {
118  paramList->validateParametersAndSetDefaults(*this->getValidParameters());
119 
120  // add in some addtional defaults that are hard to validate externally (this is because of the "disableRecursiveValidation" calls)
121 
122  if(!paramList->sublist("Initial Conditions").isType<bool>("Zero Initial Conditions"))
123  paramList->sublist("Initial Conditions").set<bool>("Zero Initial Conditions",false);
124 
125  paramList->sublist("Initial Conditions").sublist("Vector File").validateParametersAndSetDefaults(
126  getValidParameters()->sublist("Initial Conditions").sublist("Vector File"));
127 
128  this->setMyParamList(paramList);
129  }
130 
131  template<typename ScalarT>
132  Teuchos::RCP<const Teuchos::ParameterList> ModelEvaluatorFactory<ScalarT>::getValidParameters() const
133  {
134  static Teuchos::RCP<const Teuchos::ParameterList> validPL;
135  if (is_null(validPL)) {
136  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::rcp(new Teuchos::ParameterList());
137 
138  pl->sublist("Physics Blocks").disableRecursiveValidation();
139  pl->sublist("Closure Models").disableRecursiveValidation();
140  pl->sublist("Boundary Conditions").disableRecursiveValidation();
141  pl->sublist("Solution Control").disableRecursiveValidation();
142  pl->set<bool>("Use Discrete Adjoint",false);
143 
144  pl->sublist("Mesh").disableRecursiveValidation();
145 
146  pl->sublist("Initial Conditions").set<bool>("Zero Initial Conditions",false);
147  pl->sublist("Initial Conditions").sublist("Transient Parameters").disableRecursiveValidation();
148  pl->sublist("Initial Conditions").sublist("Vector File");
149  pl->sublist("Initial Conditions").sublist("Vector File").set("File Name","");
150  pl->sublist("Initial Conditions").sublist("Vector File").set<bool>("Enabled",false);
151  pl->sublist("Initial Conditions").disableRecursiveValidation();
152 
153  pl->sublist("Output").set("File Name","panzer.exo");
154  pl->sublist("Output").set("Write to Exodus",true);
155  pl->sublist("Output").sublist("Cell Average Quantities").disableRecursiveValidation();
156  pl->sublist("Output").sublist("Cell Quantities").disableRecursiveValidation();
157  pl->sublist("Output").sublist("Cell Average Vectors").disableRecursiveValidation();
158  pl->sublist("Output").sublist("Nodal Quantities").disableRecursiveValidation();
159  pl->sublist("Output").sublist("Allocate Nodal Quantities").disableRecursiveValidation();
160 
161  // Assembly sublist
162  {
163  Teuchos::ParameterList& p = pl->sublist("Assembly");
164  p.set<int>("Workset Size", 1);
165  p.set<int>("Default Integration Order",-1);
166  p.set<std::string>("Field Order","");
167  p.set<std::string>("Auxiliary Field Order","");
168  p.set<bool>("Use DOFManager FEI",false);
169  p.set<bool>("Load Balance DOFs",false);
170  p.set<bool>("Use Tpetra",false);
171  p.set<bool>("Use Epetra ME",true);
172  p.set<bool>("Lump Explicit Mass",false);
173  p.set<bool>("Constant Mass Matrix",true);
174  p.set<bool>("Apply Mass Matrix Inverse in Explicit Evaluator",true);
175  p.set<bool>("Use Conservative IMEX",false);
176  p.set<bool>("Compute Real Time Derivative",false);
177  p.set<bool>("Use Time Derivative in Explicit Model",false);
178  p.set<bool>("Compute Time Derivative at Time Step",false);
179  p.set<Teuchos::RCP<const panzer::EquationSetFactory> >("Equation Set Factory", Teuchos::null);
180  p.set<Teuchos::RCP<const panzer::ClosureModelFactory_TemplateManager<panzer::Traits> > >("Closure Model Factory", Teuchos::null);
181  p.set<Teuchos::RCP<const panzer::BCStrategyFactory> >("BC Factory",Teuchos::null);
182  p.set<std::string>("Excluded Blocks","");
183  p.sublist("ALE").disableRecursiveValidation();
184  }
185 
186  pl->sublist("Block ID to Physics ID Mapping").disableRecursiveValidation();
187  pl->sublist("Options").disableRecursiveValidation();
188  pl->sublist("Active Parameters").disableRecursiveValidation();
189  pl->sublist("Controls").disableRecursiveValidation();
190  pl->sublist("ALE").disableRecursiveValidation(); // this sucks
191  pl->sublist("User Data").disableRecursiveValidation();
192  pl->sublist("User Data").sublist("Panzer Data").disableRecursiveValidation();
193 
194  validPL = pl;
195  }
196  return validPL;
197  }
198 
199  namespace {
200  bool hasInterfaceCondition(const std::vector<panzer::BC>& bcs)
201  {
202  for (std::vector<panzer::BC>::const_iterator bcit = bcs.begin(); bcit != bcs.end(); ++bcit)
203  if (bcit->bcType() == panzer::BCT_Interface)
204  return true;
205  return false;
206  }
207 
208  Teuchos::RCP<STKConnManager>
209  getSTKConnManager(const Teuchos::RCP<panzer::ConnManager>& conn_mgr)
210  {
211  const Teuchos::RCP<STKConnManager> stk_conn_mgr =
212  Teuchos::rcp_dynamic_cast<STKConnManager>(conn_mgr);
213  TEUCHOS_TEST_FOR_EXCEPTION(stk_conn_mgr.is_null(), std::logic_error,
214  "There are interface conditions, but the connection manager"
215  " does not support the necessary connections.");
216  return stk_conn_mgr;
217  }
218 
219  void buildInterfaceConnections(const std::vector<panzer::BC>& bcs,
220  const Teuchos::RCP<panzer::ConnManager>& conn_mgr)
221  {
222  const Teuchos::RCP<STKConnManager> stk_conn_mgr = getSTKConnManager(conn_mgr);
223  for (std::vector<panzer::BC>::const_iterator bcit = bcs.begin(); bcit != bcs.end(); ++bcit)
224  if (bcit->bcType() == panzer::BCT_Interface)
225  stk_conn_mgr->associateElementsInSideset(bcit->sidesetID());
226  }
227 
228  void checkInterfaceConnections(const Teuchos::RCP<panzer::ConnManager>& conn_mgr,
229  const Teuchos::RCP<Teuchos::Comm<int> >& comm)
230  {
231  const Teuchos::RCP<STKConnManager> stk_conn_mgr = getSTKConnManager(conn_mgr);
232  std::vector<std::string> sidesets = stk_conn_mgr->checkAssociateElementsInSidesets(*comm);
233  if ( ! sidesets.empty()) {
234  std::stringstream ss;
235  ss << "Sideset IDs";
236  for (std::size_t i = 0; i < sidesets.size(); ++i)
237  ss << " " << sidesets[i];
238  ss << " did not yield associations, but these sidesets correspond to BCT_Interface BCs.";
239  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, ss.str());
240  }
241  }
242  } // namespace
243 
244  template<typename ScalarT>
245  void ModelEvaluatorFactory<ScalarT>::buildObjects(const Teuchos::RCP<const Teuchos::Comm<int> >& comm,
246  const Teuchos::RCP<panzer::GlobalData>& global_data,
247  const Teuchos::RCP<const panzer::EquationSetFactory>& eqset_factory,
248  const panzer::BCStrategyFactory & bc_factory,
250  bool meConstructionOn)
251  {
252  TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::is_null(this->getParameterList()), std::runtime_error,
253  "ParameterList must be set before objects can be built!");
254 
255  TEUCHOS_ASSERT(nonnull(comm));
256  TEUCHOS_ASSERT(nonnull(global_data));
257  TEUCHOS_ASSERT(nonnull(global_data->os));
258  TEUCHOS_ASSERT(nonnull(global_data->pl));
259 
260  // begin at the beginning...
261  m_global_data = global_data;
262 
265  // Parse input file, setup parameters
268 
269  // this function will need to be broken up eventually and probably
270  // have parts moved back into panzer. Just need to get something
271  // running.
272 
273  Teuchos::ParameterList& p = *this->getNonconstParameterList();
274 
275  // "parse" parameter list
276  Teuchos::ParameterList & mesh_params = p.sublist("Mesh");
277  Teuchos::ParameterList & assembly_params = p.sublist("Assembly");
278  Teuchos::ParameterList & solncntl_params = p.sublist("Solution Control");
279  Teuchos::ParameterList & output_list = p.sublist("Output");
280 
281  Teuchos::ParameterList & user_data_params = p.sublist("User Data");
282  Teuchos::ParameterList & panzer_data_params = user_data_params.sublist("Panzer Data");
283 
284  Teuchos::RCP<Teuchos::ParameterList> physics_block_plist = Teuchos::sublist(this->getMyNonconstParamList(),"Physics Blocks");
285 
286  // extract assembly information
287  std::size_t workset_size = Teuchos::as<std::size_t>(assembly_params.get<int>("Workset Size"));
288  std::string field_order = assembly_params.get<std::string>("Field Order"); // control nodal ordering of unknown
289  // global IDs in linear system
290  bool use_dofmanager_fei = assembly_params.get<bool>("Use DOFManager FEI"); // use FEI if true, otherwise use internal dof manager
291  bool use_load_balance = assembly_params.get<bool>("Load Balance DOFs");
292  bool useTpetra = assembly_params.get<bool>("Use Tpetra");
293  bool useThyraME = !assembly_params.get<bool>("Use Epetra ME");
294 
295  // this is weird...we are accessing the solution control to determine if things are transient
296  // it is backwards!
297  bool is_transient = solncntl_params.get<std::string>("Piro Solver") == "Rythmos" ? true : false;
298  // for pseudo-transient, we need to enable transient solver support to get time derivatives into fill
299  if (solncntl_params.get<std::string>("Piro Solver") == "NOX") {
300  if (solncntl_params.sublist("NOX").get<std::string>("Nonlinear Solver") == "Pseudo-Transient")
301  is_transient = true;
302  }
303  // for eigenvalues, we need to enable transient solver support to
304  // get time derivatives into generalized eigenvale problem
305  if (solncntl_params.get<std::string>("Piro Solver") == "LOCA") {
306  if (solncntl_params.sublist("LOCA").sublist("Stepper").get<bool>("Compute Eigenvalues"))
307  is_transient = true;
308  }
309  m_is_transient = is_transient;
310 
311  useDiscreteAdjoint = p.get<bool>("Use Discrete Adjoint");
312 
315  // Do stuff
318 
319  Teuchos::FancyOStream& fout = *global_data->os;
320 
321  // for convience cast to an MPI comm
322  const Teuchos::RCP<const Teuchos::MpiComm<int> > mpi_comm =
323  Teuchos::rcp_dynamic_cast<const Teuchos::MpiComm<int> >(comm);
324 
325  // Build mesh factory and uncommitted mesh
327 
328  Teuchos::RCP<panzer_stk::STK_MeshFactory> mesh_factory = this->buildSTKMeshFactory(mesh_params);
329  Teuchos::RCP<panzer_stk::STK_Interface> mesh = mesh_factory->buildUncommitedMesh(*(mpi_comm->getRawMpiComm()));
330  m_mesh = mesh;
331 
332  m_eqset_factory = eqset_factory;
333 
334  // setup the physcs blocks
336 
337  std::vector<Teuchos::RCP<panzer::PhysicsBlock> > physicsBlocks;
338  {
339  // setup physical mappings and boundary conditions
340  std::map<std::string,std::string> block_ids_to_physics_ids;
341  panzer::buildBlockIdToPhysicsIdMap(block_ids_to_physics_ids, p.sublist("Block ID to Physics ID Mapping"));
342 
343  // build cell ( block id -> cell topology ) mapping
344  std::map<std::string,Teuchos::RCP<const shards::CellTopology> > block_ids_to_cell_topo;
345  for(std::map<std::string,std::string>::const_iterator itr=block_ids_to_physics_ids.begin();
346  itr!=block_ids_to_physics_ids.end();++itr) {
347  block_ids_to_cell_topo[itr->first] = mesh->getCellTopology(itr->first);
348  TEUCHOS_ASSERT(block_ids_to_cell_topo[itr->first]!=Teuchos::null);
349  }
350 
351  // build physics blocks
352 
353  panzer::buildPhysicsBlocks(block_ids_to_physics_ids,
354  block_ids_to_cell_topo,
355  physics_block_plist,
356  assembly_params.get<int>("Default Integration Order"),
357  workset_size,
358  eqset_factory,
359  global_data,
360  is_transient,
361  physicsBlocks);
362  m_physics_blocks = physicsBlocks; // hold onto physics blocks for safe keeping
363  }
364 
365  // add fields automatically written through the closure model
367  addUserFieldsToMesh(*mesh,output_list);
368 
369  // finish building mesh, set required field variables and mesh bulk data
371 
372  try {
373  // this throws some exceptions, catch them as neccessary
374  this->finalizeMeshConstruction(*mesh_factory,physicsBlocks,*mpi_comm,*mesh);
375  } catch(const panzer_stk::STK_Interface::ElementBlockException & ebexp) {
376  fout << "*****************************************\n\n";
377  fout << "Element block exception, could not finalize the mesh, printing block and sideset information:\n";
378  fout.pushTab(3);
379  mesh->printMetaData(fout);
380  fout.popTab();
381  fout << std::endl;
382 
383  throw ebexp;
384  } catch(const panzer_stk::STK_Interface::SidesetException & ssexp) {
385  fout << "*****************************************\n\n";
386  fout << "Sideset exception, could not finalize the mesh, printing block and sideset information:\n";
387  fout.pushTab(3);
388  mesh->printMetaData(fout);
389  fout.popTab();
390  fout << std::endl;
391 
392  throw ssexp;
393  }
394 
395  mesh->print(fout);
396  if(p.sublist("Output").get<bool>("Write to Exodus"))
397  mesh->setupExodusFile(p.sublist("Output").get<std::string>("File Name"));
398 
399  // build a workset factory that depends on STK
401  Teuchos::RCP<panzer_stk::WorksetFactory> wkstFactory;
402  if(m_user_wkst_factory==Teuchos::null)
403  wkstFactory = Teuchos::rcp(new panzer_stk::WorksetFactory()); // build STK workset factory
404  else
405  wkstFactory = m_user_wkst_factory;
406 
407  // set workset factory mesh
408  wkstFactory->setMesh(mesh);
409 
410  // handle boundary and interface conditions
412  std::vector<panzer::BC> bcs;
413  panzer::buildBCs(bcs, p.sublist("Boundary Conditions"), global_data);
414 
415  // build the connection manager
417  Teuchos::RCP<panzer::ConnManager> conn_manager = Teuchos::rcp(new panzer_stk::STKConnManager(mesh));
418  m_conn_manager = conn_manager;
419 
420  // build DOF Manager
422 
423  Teuchos::RCP<panzer::LinearObjFactory<panzer::Traits> > linObjFactory;
424  Teuchos::RCP<panzer::GlobalIndexer> globalIndexer;
425 
426  std::string loadBalanceString = ""; // what is the load balancing information
427  bool blockedAssembly = false;
428 
429  const bool has_interface_condition = hasInterfaceCondition(bcs);
430 
431  if(panzer::BlockedDOFManagerFactory::requiresBlocking(field_order) && !useTpetra) {
432 
433  // Can't yet handle interface conditions for this system
434  TEUCHOS_TEST_FOR_EXCEPTION(has_interface_condition,
435  Teuchos::Exceptions::InvalidParameter,
436  "ERROR: Blocked Epetra systems cannot handle interface conditions.");
437 
438  // use a blocked DOF manager
439  blockedAssembly = true;
440 
441  panzer::BlockedDOFManagerFactory globalIndexerFactory;
442  globalIndexerFactory.setUseDOFManagerFEI(use_dofmanager_fei);
443 
444  Teuchos::RCP<panzer::GlobalIndexer> dofManager
445  = globalIndexerFactory.buildGlobalIndexer(mpi_comm->getRawMpiComm(),physicsBlocks,conn_manager,field_order);
446  globalIndexer = dofManager;
447 
448  Teuchos::RCP<panzer::BlockedEpetraLinearObjFactory<panzer::Traits,int> > bloLinObjFactory
450  Teuchos::rcp_dynamic_cast<panzer::BlockedDOFManager>(dofManager)));
451 
452  // parse any explicitly excluded pairs or blocks
453  const std::string excludedBlocks = assembly_params.get<std::string>("Excluded Blocks");
454  std::vector<std::string> stringPairs;
455  panzer::StringTokenizer(stringPairs,excludedBlocks,";",true);
456  for(std::size_t i=0;i<stringPairs.size();i++) {
457  std::vector<std::string> sPair;
458  std::vector<int> iPair;
459  panzer::StringTokenizer(sPair,stringPairs[i],",",true);
460  panzer::TokensToInts(iPair,sPair);
461 
462  TEUCHOS_TEST_FOR_EXCEPTION(iPair.size()!=2,std::logic_error,
463  "Input Error: The correct format for \"Excluded Blocks\" parameter in \"Assembly\" sub list is:\n"
464  " <int>,<int>; <int>,<int>; ...; <int>,<int>\n"
465  "Failure on string pair " << stringPairs[i] << "!");
466 
467  bloLinObjFactory->addExcludedPair(iPair[0],iPair[1]);
468  }
469 
470  linObjFactory = bloLinObjFactory;
471 
472  // build load balancing string for informative output
473  loadBalanceString = printUGILoadBalancingInformation(*dofManager);
474  }
475  else if(panzer::BlockedDOFManagerFactory::requiresBlocking(field_order) && useTpetra) {
476 
477  // Can't yet handle interface conditions for this system
478  TEUCHOS_TEST_FOR_EXCEPTION(has_interface_condition,
479  Teuchos::Exceptions::InvalidParameter,
480  "ERROR: Blocked Tpetra system cannot handle interface conditions.");
481 
482  // use a blocked DOF manager
483  blockedAssembly = true;
484 
485  TEUCHOS_ASSERT(!use_dofmanager_fei);
486  panzer::BlockedDOFManagerFactory globalIndexerFactory;
487  globalIndexerFactory.setUseDOFManagerFEI(false);
488 
489  Teuchos::RCP<panzer::GlobalIndexer> dofManager
490  = globalIndexerFactory.buildGlobalIndexer(mpi_comm->getRawMpiComm(),physicsBlocks,conn_manager,field_order);
491  globalIndexer = dofManager;
492 
493  Teuchos::RCP<panzer::BlockedTpetraLinearObjFactory<panzer::Traits,double,int,panzer::GlobalOrdinal> > bloLinObjFactory
495  Teuchos::rcp_dynamic_cast<panzer::BlockedDOFManager>(dofManager)));
496 
497  // parse any explicitly excluded pairs or blocks
498  const std::string excludedBlocks = assembly_params.get<std::string>("Excluded Blocks");
499  std::vector<std::string> stringPairs;
500  panzer::StringTokenizer(stringPairs,excludedBlocks,";",true);
501  for(std::size_t i=0;i<stringPairs.size();i++) {
502  std::vector<std::string> sPair;
503  std::vector<int> iPair;
504  panzer::StringTokenizer(sPair,stringPairs[i],",",true);
505  panzer::TokensToInts(iPair,sPair);
506 
507  TEUCHOS_TEST_FOR_EXCEPTION(iPair.size()!=2,std::logic_error,
508  "Input Error: The correct format for \"Excluded Blocks\" parameter in \"Assembly\" sub list is:\n"
509  " <int>,<int>; <int>,<int>; ...; <int>,<int>\n"
510  "Failure on string pair " << stringPairs[i] << "!");
511 
512  bloLinObjFactory->addExcludedPair(iPair[0],iPair[1]);
513  }
514 
515  linObjFactory = bloLinObjFactory;
516 
517  // build load balancing string for informative output
518  loadBalanceString = printUGILoadBalancingInformation(*dofManager);
519  }
520  else if(useTpetra) {
521 
522  if (has_interface_condition)
523  buildInterfaceConnections(bcs, conn_manager);
524 
525  // use a flat DOF manager
526 
527  TEUCHOS_ASSERT(!use_dofmanager_fei);
528  panzer::DOFManagerFactory globalIndexerFactory;
529  globalIndexerFactory.setUseDOFManagerFEI(false);
530  globalIndexerFactory.setUseTieBreak(use_load_balance);
531  Teuchos::RCP<panzer::GlobalIndexer> dofManager
532  = globalIndexerFactory.buildGlobalIndexer(mpi_comm->getRawMpiComm(),physicsBlocks,conn_manager,field_order);
533  globalIndexer = dofManager;
534 
535  if (has_interface_condition)
536  checkInterfaceConnections(conn_manager, dofManager->getComm());
537 
538  TEUCHOS_ASSERT(!useDiscreteAdjoint); // safety check
539  linObjFactory = Teuchos::rcp(new panzer::TpetraLinearObjFactory<panzer::Traits,double,int,panzer::GlobalOrdinal>(mpi_comm,dofManager));
540 
541  // build load balancing string for informative output
542  loadBalanceString = printUGILoadBalancingInformation(*dofManager);
543  }
544  else {
545 
546  if (has_interface_condition)
547  buildInterfaceConnections(bcs, conn_manager);
548 
549  // use a flat DOF manager
550  panzer::DOFManagerFactory globalIndexerFactory;
551  globalIndexerFactory.setUseDOFManagerFEI(use_dofmanager_fei);
552  globalIndexerFactory.setUseTieBreak(use_load_balance);
553  globalIndexerFactory.setUseNeighbors(has_interface_condition);
554  Teuchos::RCP<panzer::GlobalIndexer> dofManager
555  = globalIndexerFactory.buildGlobalIndexer(mpi_comm->getRawMpiComm(),physicsBlocks,conn_manager,
556  field_order);
557  globalIndexer = dofManager;
558 
559  if (has_interface_condition)
560  checkInterfaceConnections(conn_manager, dofManager->getComm());
561 
562  linObjFactory = Teuchos::rcp(new panzer::BlockedEpetraLinearObjFactory<panzer::Traits,int>(mpi_comm,dofManager,useDiscreteAdjoint));
563 
564  // build load balancing string for informative output
565  loadBalanceString = printUGILoadBalancingInformation(*dofManager);
566  }
567 
568  TEUCHOS_ASSERT(globalIndexer!=Teuchos::null);
569  TEUCHOS_ASSERT(linObjFactory!=Teuchos::null);
570  m_global_indexer = globalIndexer;
571  m_lin_obj_factory = linObjFactory;
572  m_blockedAssembly = blockedAssembly;
573 
574  // print out load balancing information
575  fout << "Degree of freedom load balancing: " << loadBalanceString << std::endl;
576 
577  // build worksets
579 
580  // build up needs array for workset container
581  std::map<std::string,panzer::WorksetNeeds> needs;
582  for(std::size_t i=0;i<physicsBlocks.size();i++)
583  needs[physicsBlocks[i]->elementBlockID()] = physicsBlocks[i]->getWorksetNeeds();
584 
585  Teuchos::RCP<panzer::WorksetContainer> wkstContainer // attach it to a workset container (uses lazy evaluation)
586  = Teuchos::rcp(new panzer::WorksetContainer(wkstFactory,needs));
587 
588  wkstContainer->setWorksetSize(workset_size);
589  wkstContainer->setGlobalIndexer(globalIndexer); // set the global indexer so the orientations are evaluated
590 
591  m_wkstContainer = wkstContainer;
592 
593  // find max number of worksets
594  std::size_t max_wksets = 0;
595  for(std::size_t pb=0;pb<physicsBlocks.size();pb++) {
596  const panzer::WorksetDescriptor wd = panzer::blockDescriptor(physicsBlocks[pb]->elementBlockID());
597  Teuchos::RCP< std::vector<panzer::Workset> >works = wkstContainer->getWorksets(wd);
598  max_wksets = std::max(max_wksets,works->size());
599  }
600  user_data_params.set<std::size_t>("Max Worksets",max_wksets);
601  wkstContainer->clear();
602 
603  // Setup lagrangian type coordinates
605 
606  // see if field coordinates are required, if so reset the workset container
607  // and set the coordinates to be associated with a field in the mesh
608  useDynamicCoordinates_ = false;
609  for(std::size_t pb=0;pb<physicsBlocks.size();pb++) {
610  if(physicsBlocks[pb]->getCoordinateDOFs().size()>0) {
611  mesh->setUseFieldCoordinates(true);
612  useDynamicCoordinates_ = true;
613  wkstContainer->clear(); // this serves to refresh the worksets
614  // and put in new coordinates
615  break;
616  }
617  }
618 
619  // Add mesh objects to user data to make available to user ctors
621 
622  panzer_data_params.set("STK Mesh", mesh);
623  panzer_data_params.set("DOF Manager", globalIndexer);
624  panzer_data_params.set("Linear Object Factory", linObjFactory);
625 
626  // If user requested it, short circuit model construction
628 
629  if(!meConstructionOn)
630  return;
631 
632  // Setup active parameters
634 
635  std::vector<Teuchos::RCP<Teuchos::Array<std::string> > > p_names;
636  std::vector<Teuchos::RCP<Teuchos::Array<double> > > p_values;
637  if (p.isSublist("Active Parameters")) {
638  Teuchos::ParameterList& active_params = p.sublist("Active Parameters");
639 
640  int num_param_vecs = active_params.get<int>("Number of Parameter Vectors",0);
641  p_names.resize(num_param_vecs);
642  p_values.resize(num_param_vecs);
643  for (int i=0; i<num_param_vecs; i++) {
644  std::stringstream ss;
645  ss << "Parameter Vector " << i;
646  Teuchos::ParameterList& pList = active_params.sublist(ss.str());
647  int numParameters = pList.get<int>("Number");
648  TEUCHOS_TEST_FOR_EXCEPTION(numParameters == 0,
649  Teuchos::Exceptions::InvalidParameter,
650  std::endl << "Error! panzer::ModelEvaluator::ModelEvaluator(): " <<
651  "Parameter vector " << i << " has zero parameters!" << std::endl);
652  p_names[i] =
653  Teuchos::rcp(new Teuchos::Array<std::string>(numParameters));
654  p_values[i] =
655  Teuchos::rcp(new Teuchos::Array<double>(numParameters));
656  for (int j=0; j<numParameters; j++) {
657  std::stringstream ss2;
658  ss2 << "Parameter " << j;
659  (*p_names[i])[j] = pList.get<std::string>(ss2.str());
660  ss2.str("");
661 
662  ss2 << "Initial Value " << j;
663  (*p_values[i])[j] = pList.get<double>(ss2.str());
664 
665  // this is a band-aid/hack to make sure parameters are registered before they are accessed
666  panzer::registerScalarParameter((*p_names[i])[j],*global_data->pl,(*p_values[i])[j]);
667  }
668  }
669  }
670 
671  // setup the closure model for automatic writing (during residual/jacobian update)
673 
674  panzer_stk::IOClosureModelFactory_TemplateBuilder<panzer::Traits> io_cm_builder(user_cm_factory,mesh,output_list);
676  cm_factory.buildObjects(io_cm_builder);
677 
678  // setup field manager build
680 
681  Teuchos::RCP<panzer::FieldManagerBuilder> fmb;
682  {
683  bool write_dot_files = p.sublist("Options").get("Write Volume Assembly Graphs",false);
684  std::string dot_file_prefix = p.sublist("Options").get("Volume Assembly Graph Prefix","Panzer_AssemblyGraph");
685  bool write_fm_files = p.sublist("Options").get("Write Field Manager Files",false);
686  std::string fm_file_prefix = p.sublist("Options").get("Field Manager File Prefix","Panzer_AssemblyGraph");
687 
688  // Allow users to override inputs via runtime env
689  {
690  auto check_write_dag = std::getenv("PANZER_WRITE_DAG");
691  if (check_write_dag != nullptr) {
692  write_dot_files = true;
693  write_fm_files = true;
694  }
695  }
696 
697  fmb = buildFieldManagerBuilder(wkstContainer,physicsBlocks,bcs,*eqset_factory,bc_factory,cm_factory,
698  user_cm_factory,p.sublist("Closure Models"),*linObjFactory,user_data_params,
699  write_dot_files,dot_file_prefix,
700  write_fm_files,fm_file_prefix);
701  }
702 
703  // build response library
705 
706  m_response_library = Teuchos::rcp(new panzer::ResponseLibrary<panzer::Traits>(wkstContainer,globalIndexer,linObjFactory));
707 
708  {
709  bool write_dot_files = false;
710  std::string prefix = "Panzer_ResponseGraph_";
711  write_dot_files = p.sublist("Options").get("Write Volume Response Graphs",write_dot_files);
712  prefix = p.sublist("Options").get("Volume Response Graph Prefix",prefix);
713 
714  Teuchos::ParameterList user_data(p.sublist("User Data"));
715  user_data.set<int>("Workset Size",workset_size);
716  }
717 
718  // Setup solver factory
720 
721  Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<double> > lowsFactory =
722  buildLOWSFactory(blockedAssembly,globalIndexer,conn_manager,mesh,mpi_comm);
723 
724  // Setup physics model evaluator
726 
727  double t_init = 0.0;
728  if(is_transient)
729  t_init = this->getInitialTime(p.sublist("Initial Conditions").sublist("Transient Parameters"), *mesh);
730 
731  if(blockedAssembly || useTpetra) // override the user request
732  useThyraME = true;
733 
734  Teuchos::RCP<Thyra::ModelEvaluatorDefaultBase<double> > thyra_me
735  = buildPhysicsModelEvaluator(useThyraME, // blockedAssembly || useTpetra, // this determines if a Thyra or Epetra ME is used
736  fmb,
737  m_response_library,
738  linObjFactory,
739  p_names,
740  p_values,
741  lowsFactory,
742  global_data,
743  is_transient,
744  t_init);
745 
746  // Setup initial conditions
748 
749  {
750  // Create closure model list for use in defining initial conditions
751  // For now just remove Global MMS Parameters, if it exists
752  const Teuchos::ParameterList& models = p.sublist("Closure Models");
753  Teuchos::ParameterList cl_models(models.name());
754  for (Teuchos::ParameterList::ConstIterator model_it=models.begin();
755  model_it!=models.end(); ++model_it) {
756  std::string key = model_it->first;
757  if (model_it->first != "Global MMS Parameters")
758  cl_models.setEntry(key,model_it->second);
759  }
760  bool write_dot_files = false;
761  std::string prefix = "Panzer_AssemblyGraph_";
762  setupInitialConditions(*thyra_me,*wkstContainer,physicsBlocks,user_cm_factory,*linObjFactory,
763  cl_models,
764  p.sublist("Initial Conditions"),
765  p.sublist("User Data"),
766  p.sublist("Options").get("Write Volume Assembly Graphs",write_dot_files),
767  p.sublist("Options").get("Volume Assembly Graph Prefix",prefix));
768  }
769 
770  // Write the IC vector into the STK mesh: use response library
772  writeInitialConditions(*thyra_me,physicsBlocks,wkstContainer,globalIndexer,linObjFactory,mesh,user_cm_factory,
773  p.sublist("Closure Models"),
774  p.sublist("User Data"),workset_size);
775 
776  m_physics_me = thyra_me;
777  }
778 
779  template<typename ScalarT>
781  addUserFieldsToMesh(panzer_stk::STK_Interface & mesh,const Teuchos::ParameterList & output_list) const
782  {
783  // register cell averaged scalar fields
784  const Teuchos::ParameterList & cellAvgQuants = output_list.sublist("Cell Average Quantities");
785  for(Teuchos::ParameterList::ConstIterator itr=cellAvgQuants.begin();
786  itr!=cellAvgQuants.end();++itr) {
787  const std::string & blockId = itr->first;
788  const std::string & fields = Teuchos::any_cast<std::string>(itr->second.getAny());
789  std::vector<std::string> tokens;
790 
791  // break up comma seperated fields
792  panzer::StringTokenizer(tokens,fields,",",true);
793 
794  for(std::size_t i=0;i<tokens.size();i++)
795  mesh.addCellField(tokens[i],blockId);
796  }
797 
798  // register cell averaged components of vector fields
799  // just allocate space for the fields here. The actual calculation and writing
800  // are done by panzer_stk::ScatterCellAvgVector.
801  const Teuchos::ParameterList & cellAvgVectors = output_list.sublist("Cell Average Vectors");
802  for(Teuchos::ParameterList::ConstIterator itr = cellAvgVectors.begin();
803  itr != cellAvgVectors.end(); ++itr) {
804  const std::string & blockId = itr->first;
805  const std::string & fields = Teuchos::any_cast<std::string>(itr->second.getAny());
806  std::vector<std::string> tokens;
807 
808  // break up comma seperated fields
809  panzer::StringTokenizer(tokens,fields,",",true);
810 
811  for(std::size_t i = 0; i < tokens.size(); i++) {
812  std::string d_mod[3] = {"X","Y","Z"};
813  for(std::size_t d = 0; d < mesh.getDimension(); d++)
814  mesh.addCellField(tokens[i]+d_mod[d],blockId);
815  }
816  }
817 
818  // register cell quantities
819  const Teuchos::ParameterList & cellQuants = output_list.sublist("Cell Quantities");
820  for(Teuchos::ParameterList::ConstIterator itr=cellQuants.begin();
821  itr!=cellQuants.end();++itr) {
822  const std::string & blockId = itr->first;
823  const std::string & fields = Teuchos::any_cast<std::string>(itr->second.getAny());
824  std::vector<std::string> tokens;
825 
826  // break up comma seperated fields
827  panzer::StringTokenizer(tokens,fields,",",true);
828 
829  for(std::size_t i=0;i<tokens.size();i++)
830  mesh.addCellField(tokens[i],blockId);
831  }
832 
833  // register ndoal quantities
834  const Teuchos::ParameterList & nodalQuants = output_list.sublist("Nodal Quantities");
835  for(Teuchos::ParameterList::ConstIterator itr=nodalQuants.begin();
836  itr!=nodalQuants.end();++itr) {
837  const std::string & blockId = itr->first;
838  const std::string & fields = Teuchos::any_cast<std::string>(itr->second.getAny());
839  std::vector<std::string> tokens;
840 
841  // break up comma seperated fields
842  panzer::StringTokenizer(tokens,fields,",",true);
843 
844  for(std::size_t i=0;i<tokens.size();i++)
845  mesh.addSolutionField(tokens[i],blockId);
846  }
847 
848  const Teuchos::ParameterList & allocNodalQuants = output_list.sublist("Allocate Nodal Quantities");
849  for(Teuchos::ParameterList::ConstIterator itr=allocNodalQuants.begin();
850  itr!=allocNodalQuants.end();++itr) {
851  const std::string & blockId = itr->first;
852  const std::string & fields = Teuchos::any_cast<std::string>(itr->second.getAny());
853  std::vector<std::string> tokens;
854 
855  // break up comma seperated fields
856  panzer::StringTokenizer(tokens,fields,",",true);
857 
858  for(std::size_t i=0;i<tokens.size();i++)
859  mesh.addSolutionField(tokens[i],blockId);
860  }
861  }
862 
863  template<typename ScalarT>
866  panzer::WorksetContainer & wkstContainer,
867  const std::vector<Teuchos::RCP<panzer::PhysicsBlock> >& physicsBlocks,
870  const Teuchos::ParameterList & closure_pl,
871  const Teuchos::ParameterList & initial_cond_pl,
872  const Teuchos::ParameterList & user_data_pl,
873  bool write_dot_files,const std::string & dot_file_prefix) const
874  {
875  using Teuchos::RCP;
876 
877  Thyra::ModelEvaluatorBase::InArgs<double> nomValues = model.getNominalValues();
878  Teuchos::RCP<Thyra::VectorBase<double> > x_vec = Teuchos::rcp_const_cast<Thyra::VectorBase<double> >(nomValues.get_x());
879 
880  if(initial_cond_pl.get<bool>("Zero Initial Conditions")) {
881  // zero out the x vector
882  Thyra::assign(x_vec.ptr(),0.0);
883  }
884  else if(!initial_cond_pl.sublist("Vector File").get<bool>("Enabled")) {
885  // read from exodus, or compute using field managers
886 
887  std::map<std::string, Teuchos::RCP< PHX::FieldManager<panzer::Traits> > > phx_ic_field_managers;
888 
890  physicsBlocks,
891  cm_factory,
892  closure_pl,
893  initial_cond_pl,
894  lof,
895  user_data_pl,
896  write_dot_files,
897  dot_file_prefix,
898  phx_ic_field_managers);
899 /*
900  panzer::setupInitialConditionFieldManagers(wkstContainer,
901  physicsBlocks,
902  cm_factory,
903  initial_cond_pl,
904  lof,
905  user_data_pl,
906  write_dot_files,
907  dot_file_prefix,
908  phx_ic_field_managers);
909 */
910 
911  // set the vector to be filled
912  Teuchos::RCP<panzer::LinearObjContainer> loc = lof.buildLinearObjContainer();
913  Teuchos::RCP<panzer::ThyraObjContainer<double> > tloc = Teuchos::rcp_dynamic_cast<panzer::ThyraObjContainer<double> >(loc);
914  tloc->set_x_th(x_vec);
915 
916  panzer::evaluateInitialCondition(wkstContainer, phx_ic_field_managers, loc, lof, 0.0);
917  }
918  else {
919  const std::string & vectorFile = initial_cond_pl.sublist("Vector File").get<std::string>("File Name");
920  TEUCHOS_TEST_FOR_EXCEPTION(vectorFile=="",std::runtime_error,
921  "If \"Read From Vector File\" is true, then parameter \"Vector File\" cannot be the empty string.");
922 
923  // set the vector to be filled
924  Teuchos::RCP<panzer::LinearObjContainer> loc = lof.buildLinearObjContainer();
925  Teuchos::RCP<panzer::ThyraObjContainer<double> > tloc = Teuchos::rcp_dynamic_cast<panzer::ThyraObjContainer<double> >(loc);
926  tloc->set_x_th(x_vec);
927 
928  // read the vector
929  lof.readVector(vectorFile,*loc,panzer::LinearObjContainer::X);
930  }
931  }
932 
933  template<typename ScalarT>
936  const std::vector<Teuchos::RCP<panzer::PhysicsBlock> >& physicsBlocks,
937  const Teuchos::RCP<panzer::WorksetContainer> & wc,
938  const Teuchos::RCP<const panzer::GlobalIndexer> & ugi,
939  const Teuchos::RCP<const panzer::LinearObjFactory<panzer::Traits> > & lof,
940  const Teuchos::RCP<panzer_stk::STK_Interface> & mesh,
942  const Teuchos::ParameterList & closure_model_pl,
943  const Teuchos::ParameterList & user_data_pl,
944  int workset_size) const
945  {
946  Teuchos::RCP<panzer::LinearObjContainer> loc = lof->buildLinearObjContainer();
947  Teuchos::RCP<panzer::ThyraObjContainer<double> > tloc = Teuchos::rcp_dynamic_cast<panzer::ThyraObjContainer<double> >(loc);
948  tloc->set_x_th(Teuchos::rcp_const_cast<Thyra::VectorBase<double> >(model.getNominalValues().get_x()));
949 
950  Teuchos::RCP<panzer::ResponseLibrary<panzer::Traits> > solnWriter
951  = initializeSolnWriterResponseLibrary(wc,ugi,lof,mesh);
952 
953  {
954  Teuchos::ParameterList user_data(user_data_pl);
955  user_data.set<int>("Workset Size",workset_size);
956 
957  finalizeSolnWriterResponseLibrary(*solnWriter,physicsBlocks,cm_factory,closure_model_pl,workset_size,user_data);
958  }
959 
960  // initialize the assembly container
962  ae_inargs.container_ = loc;
963  ae_inargs.ghostedContainer_ = lof->buildGhostedLinearObjContainer();
964  ae_inargs.alpha = 0.0;
965  ae_inargs.beta = 1.0;
966  ae_inargs.evaluate_transient_terms = false;
967 
968  // initialize the ghosted container
969  lof->initializeGhostedContainer(panzer::LinearObjContainer::X,*ae_inargs.ghostedContainer_);
970 
971  // do import
972  lof->globalToGhostContainer(*ae_inargs.container_,*ae_inargs.ghostedContainer_,panzer::LinearObjContainer::X);
973 
974  // fill STK mesh objects
975  solnWriter->addResponsesToInArgs<panzer::Traits::Residual>(ae_inargs);
976  solnWriter->evaluate<panzer::Traits::Residual>(ae_inargs);
977  }
978 
980  template<typename ScalarT>
981  Teuchos::RCP<panzer_stk::STK_MeshFactory> ModelEvaluatorFactory<ScalarT>::buildSTKMeshFactory(const Teuchos::ParameterList & mesh_params) const
982  {
983  Teuchos::RCP<panzer_stk::STK_MeshFactory> mesh_factory;
984 
985  // first contruct the mesh factory
986  if (mesh_params.get<std::string>("Source") == "Exodus File") {
987  mesh_factory = Teuchos::rcp(new panzer_stk::STK_ExodusReaderFactory());
988  mesh_factory->setParameterList(Teuchos::rcp(new Teuchos::ParameterList(mesh_params.sublist("Exodus File"))));
989  }
990  else if (mesh_params.get<std::string>("Source") == "Pamgen Mesh") {
991  mesh_factory = Teuchos::rcp(new panzer_stk::STK_ExodusReaderFactory());
992  Teuchos::RCP<Teuchos::ParameterList> pamgenList = Teuchos::rcp(new Teuchos::ParameterList(mesh_params.sublist("Pamgen Mesh")));
993  pamgenList->set("File Type","Pamgen"); // For backwards compatibility when pamgen had separate factory from exodus
994  mesh_factory->setParameterList(pamgenList);
995  }
996  else if (mesh_params.get<std::string>("Source") == "Inline Mesh") {
997 
998  int dimension = mesh_params.sublist("Inline Mesh").get<int>("Mesh Dimension");
999  std::string typeStr = "";
1000  if(mesh_params.sublist("Inline Mesh").isParameter("Type"))
1001  typeStr = mesh_params.sublist("Inline Mesh").get<std::string>("Type");
1002 
1003  if (dimension == 1) {
1004  mesh_factory = Teuchos::rcp(new panzer_stk::LineMeshFactory);
1005  Teuchos::RCP<Teuchos::ParameterList> in_mesh = Teuchos::rcp(new Teuchos::ParameterList);
1006  *in_mesh = mesh_params.sublist("Inline Mesh").sublist("Mesh Factory Parameter List");
1007  mesh_factory->setParameterList(in_mesh);
1008  }
1009  else if (dimension == 2 && typeStr=="Tri") {
1010  mesh_factory = Teuchos::rcp(new panzer_stk::SquareTriMeshFactory);
1011  Teuchos::RCP<Teuchos::ParameterList> in_mesh = Teuchos::rcp(new Teuchos::ParameterList);
1012  *in_mesh = mesh_params.sublist("Inline Mesh").sublist("Mesh Factory Parameter List");
1013  mesh_factory->setParameterList(in_mesh);
1014  }
1015  else if (dimension == 2) {
1016  mesh_factory = Teuchos::rcp(new panzer_stk::SquareQuadMeshFactory);
1017  Teuchos::RCP<Teuchos::ParameterList> in_mesh = Teuchos::rcp(new Teuchos::ParameterList);
1018  *in_mesh = mesh_params.sublist("Inline Mesh").sublist("Mesh Factory Parameter List");
1019  mesh_factory->setParameterList(in_mesh);
1020  }
1021  else if (dimension == 3 && typeStr=="Tet") {
1022  mesh_factory = Teuchos::rcp(new panzer_stk::CubeTetMeshFactory);
1023  Teuchos::RCP<Teuchos::ParameterList> in_mesh = Teuchos::rcp(new Teuchos::ParameterList);
1024  *in_mesh = mesh_params.sublist("Inline Mesh").sublist("Mesh Factory Parameter List");
1025  mesh_factory->setParameterList(in_mesh);
1026  }
1027  else if(dimension == 3) {
1028  mesh_factory = Teuchos::rcp(new panzer_stk::CubeHexMeshFactory);
1029  Teuchos::RCP<Teuchos::ParameterList> in_mesh = Teuchos::rcp(new Teuchos::ParameterList);
1030  *in_mesh = mesh_params.sublist("Inline Mesh").sublist("Mesh Factory Parameter List");
1031  mesh_factory->setParameterList(in_mesh);
1032  }
1033  else if(dimension==4) { // not really "dimension==4" simply a flag to try this other mesh for testing
1034  mesh_factory = Teuchos::rcp(new panzer_stk::MultiBlockMeshFactory);
1035  Teuchos::RCP<Teuchos::ParameterList> in_mesh = Teuchos::rcp(new Teuchos::ParameterList);
1036  *in_mesh = mesh_params.sublist("Inline Mesh").sublist("Mesh Factory Parameter List");
1037  mesh_factory->setParameterList(in_mesh);
1038  }
1039  }
1040  else if (mesh_params.get<std::string>("Source") == "Custom Mesh") {
1041  mesh_factory = Teuchos::rcp(new panzer_stk::CustomMeshFactory());
1042  mesh_factory->setParameterList(Teuchos::rcp(new Teuchos::ParameterList(mesh_params.sublist("Custom Mesh"))));
1043  }
1044  else {
1045  // throw a runtime exception for invalid parameter values
1046  }
1047 
1048 
1049  // get rebalancing parameters
1050  if(mesh_params.isSublist("Rebalance")) {
1051  const Teuchos::ParameterList & rebalance = mesh_params.sublist("Rebalance");
1052 
1053  // check to see if its enabled
1054  bool enabled = false;
1055  if(rebalance.isType<bool>("Enabled"))
1056  enabled = rebalance.get<bool>("Enabled");
1057 
1058  // we can also use a list description of what to load balance
1059  Teuchos::RCP<Teuchos::ParameterList> rebalanceCycles;
1060  if(enabled && rebalance.isSublist("Cycles"))
1061  rebalanceCycles = Teuchos::rcp(new Teuchos::ParameterList(rebalance.sublist("Cycles")));
1062 
1063  // setup rebalancing as neccessary
1064  mesh_factory->enableRebalance(enabled,rebalanceCycles);
1065  }
1066 
1067  return mesh_factory;
1068  }
1069 
1070  template<typename ScalarT>
1072  const std::vector<Teuchos::RCP<panzer::PhysicsBlock> > & physicsBlocks,
1073  const Teuchos::MpiComm<int> mpi_comm,
1074  STK_Interface & mesh) const
1075  {
1076  // finish building mesh, set required field variables and mesh bulk data
1077  {
1078  std::vector<Teuchos::RCP<panzer::PhysicsBlock> >::const_iterator physIter;
1079  for(physIter=physicsBlocks.begin();physIter!=physicsBlocks.end();++physIter) {
1080  // what is the block weight for this element block?
1081  double blockWeight = 0.0;
1082 
1083  Teuchos::RCP<const panzer::PhysicsBlock> pb = *physIter;
1084  const std::vector<panzer::StrPureBasisPair> & blockFields = pb->getProvidedDOFs();
1085  const std::vector<std::vector<std::string> > & coordinateDOFs = pb->getCoordinateDOFs();
1086  // these are treated specially
1087 
1088  // insert all fields into a set
1089  std::set<panzer::StrPureBasisPair,panzer::StrPureBasisComp> fieldNames;
1090  fieldNames.insert(blockFields.begin(),blockFields.end());
1091 
1092  // Now we will set up the coordinate fields (make sure to remove
1093  // the DOF fields)
1094  {
1095  std::set<std::string> fields_to_remove;
1096 
1097  // add mesh coordinate fields, setup their removal from fieldNames
1098  // set to prevent duplication
1099  for(std::size_t i=0;i<coordinateDOFs.size();i++) {
1100  mesh.addMeshCoordFields(pb->elementBlockID(),coordinateDOFs[i],"DISPL");
1101  for(std::size_t j=0;j<coordinateDOFs[i].size();j++)
1102  fields_to_remove.insert(coordinateDOFs[i][j]);
1103  }
1104 
1105  // remove the already added coordinate fields
1106  std::set<std::string>::const_iterator rmItr;
1107  for (rmItr=fields_to_remove.begin();rmItr!=fields_to_remove.end();++rmItr)
1108  fieldNames.erase(fieldNames.find(panzer::StrPureBasisPair(*rmItr,Teuchos::null)));
1109  }
1110 
1111  // add basis to DOF manager: block specific
1112  std::set<panzer::StrPureBasisPair,panzer::StrPureBasisComp>::const_iterator fieldItr;
1113  for (fieldItr=fieldNames.begin();fieldItr!=fieldNames.end();++fieldItr) {
1114 
1115  if(fieldItr->second->isScalarBasis() &&
1116  fieldItr->second->getElementSpace()==panzer::PureBasis::CONST) {
1117  mesh.addCellField(fieldItr->first,pb->elementBlockID());
1118  }
1119  else if(fieldItr->second->isScalarBasis()) {
1120  mesh.addSolutionField(fieldItr->first,pb->elementBlockID());
1121  }
1122  else if(fieldItr->second->isVectorBasis()) {
1123  std::string d_mod[3] = {"X","Y","Z"};
1124  for(int d=0;d<fieldItr->second->dimension();d++)
1125  mesh.addCellField(fieldItr->first+d_mod[d],pb->elementBlockID());
1126  }
1127  else { TEUCHOS_ASSERT(false); }
1128 
1129  blockWeight += double(fieldItr->second->cardinality());
1130  }
1131 
1132  // set the compute block weight (this is the sum of the cardinality of all basis
1133  // functions on this block
1134  mesh.setBlockWeight(pb->elementBlockID(),blockWeight);
1135  }
1136 
1137  mesh_factory.completeMeshConstruction(mesh,*(mpi_comm.getRawMpiComm()));
1138  }
1139  }
1140 
1141 
1142  template<typename ScalarT>
1143  Teuchos::RCP<Thyra::ModelEvaluator<ScalarT> > ModelEvaluatorFactory<ScalarT>::getPhysicsModelEvaluator()
1144  {
1145  TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::is_null(m_physics_me), std::runtime_error,
1146  "Objects are not built yet! Please call buildObjects() member function.");
1147  return m_physics_me;
1148  }
1149 
1150  template<typename ScalarT>
1151  void ModelEvaluatorFactory<ScalarT>::setNOXObserverFactory(const Teuchos::RCP<const panzer_stk::NOXObserverFactory>& nox_observer_factory)
1152  {
1153  m_nox_observer_factory = nox_observer_factory;
1154  }
1155 
1156  template<typename ScalarT>
1157  void ModelEvaluatorFactory<ScalarT>::setRythmosObserverFactory(const Teuchos::RCP<const panzer_stk::RythmosObserverFactory>& rythmos_observer_factory)
1158  {
1159  m_rythmos_observer_factory = rythmos_observer_factory;
1160  }
1161 
1162  template<typename ScalarT>
1163  void ModelEvaluatorFactory<ScalarT>::setUserWorksetFactory(Teuchos::RCP<panzer_stk::WorksetFactory>& user_wkst_factory)
1164  {
1165  m_user_wkst_factory = user_wkst_factory;
1166  }
1167 
1168  template<typename ScalarT>
1169  Teuchos::RCP<Thyra::ModelEvaluator<ScalarT> > ModelEvaluatorFactory<ScalarT>::getResponseOnlyModelEvaluator()
1170  {
1171  if(m_rome_me==Teuchos::null)
1172  m_rome_me = buildResponseOnlyModelEvaluator(m_physics_me,m_global_data);
1173 
1174  return m_rome_me;
1175  }
1176 
1177  template<typename ScalarT>
1178  Teuchos::RCP<Thyra::ModelEvaluator<ScalarT> > ModelEvaluatorFactory<ScalarT>::
1180  const Teuchos::RCP<panzer::GlobalData>& global_data,
1181  const Teuchos::RCP<Piro::RythmosSolver<ScalarT> > rythmosSolver,
1182  const Teuchos::Ptr<const panzer_stk::NOXObserverFactory> & in_nox_observer_factory,
1183  const Teuchos::Ptr<const panzer_stk::RythmosObserverFactory> & in_rythmos_observer_factory
1184  )
1185  {
1186  using Teuchos::is_null;
1187  using Teuchos::Ptr;
1188 
1189  TEUCHOS_TEST_FOR_EXCEPTION(is_null(m_lin_obj_factory), std::runtime_error,
1190  "Objects are not built yet! Please call buildObjects() member function.");
1191  TEUCHOS_TEST_FOR_EXCEPTION(is_null(m_global_indexer), std::runtime_error,
1192  "Objects are not built yet! Please call buildObjects() member function.");
1193  TEUCHOS_TEST_FOR_EXCEPTION(is_null(m_mesh), std::runtime_error,
1194  "Objects are not built yet! Please call buildObjects() member function.");
1195  Teuchos::Ptr<const panzer_stk::NOXObserverFactory> nox_observer_factory
1196  = is_null(in_nox_observer_factory) ? m_nox_observer_factory.ptr() : in_nox_observer_factory;
1197  Teuchos::Ptr<const panzer_stk::RythmosObserverFactory> rythmos_observer_factory
1198  = is_null(in_rythmos_observer_factory) ? m_rythmos_observer_factory.ptr() : in_rythmos_observer_factory;
1199 
1200  Teuchos::ParameterList& p = *this->getNonconstParameterList();
1201  Teuchos::ParameterList & solncntl_params = p.sublist("Solution Control");
1202  Teuchos::RCP<Teuchos::ParameterList> piro_params = Teuchos::rcp(new Teuchos::ParameterList(solncntl_params));
1203  Teuchos::RCP<Thyra::ModelEvaluatorDefaultBase<double> > piro;
1204 
1205  std::string solver = solncntl_params.get<std::string>("Piro Solver");
1206  Teuchos::RCP<Thyra::ModelEvaluatorDefaultBase<double> > thyra_me_db
1207  = Teuchos::rcp_dynamic_cast<Thyra::ModelEvaluatorDefaultBase<double> >(thyra_me);
1208  if ( (solver=="NOX") || (solver == "LOCA") ) {
1209 
1210  TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::is_null(nox_observer_factory), std::runtime_error,
1211  "No NOX obersver built! Please call setNOXObserverFactory() member function if you plan to use a NOX solver.");
1212 
1213  Teuchos::RCP<NOX::Abstract::PrePostOperator> ppo = nox_observer_factory->buildNOXObserver(m_mesh,m_global_indexer,m_lin_obj_factory);
1214  piro_params->sublist("NOX").sublist("Solver Options").set("User Defined Pre/Post Operator", ppo);
1215 
1216  if (solver=="NOX")
1217  piro = Teuchos::rcp(new Piro::NOXSolver<double>(piro_params,
1218  Teuchos::rcp_dynamic_cast<Thyra::ModelEvaluatorDefaultBase<double> >(thyra_me_db)));
1219  else if (solver == "LOCA")
1220  piro = Teuchos::rcp(new Piro::LOCASolver<double>(piro_params,
1221  Teuchos::rcp_dynamic_cast<Thyra::ModelEvaluatorDefaultBase<double> >(thyra_me_db),
1222  Teuchos::null));
1223  TEUCHOS_ASSERT(nonnull(piro));
1224 
1225  // override printing to use panzer ostream
1226  piro_params->sublist("NOX").sublist("Printing").set<Teuchos::RCP<std::ostream> >("Output Stream",global_data->os);
1227  piro_params->sublist("NOX").sublist("Printing").set<Teuchos::RCP<std::ostream> >("Error Stream",global_data->os);
1228  piro_params->sublist("NOX").sublist("Printing").set<int>("Output Processor",global_data->os->getOutputToRootOnly());
1229  }
1230  else if (solver=="Rythmos") {
1231 
1232  TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::is_null(rythmos_observer_factory), std::runtime_error,
1233  "No NOX obersver built! Please call setrythmosObserverFactory() member function if you plan to use a Rythmos solver.");
1234 
1235  // install the nox observer
1236  if(rythmos_observer_factory->useNOXObserver()) {
1237  Teuchos::RCP<NOX::Abstract::PrePostOperator> ppo = nox_observer_factory->buildNOXObserver(m_mesh,m_global_indexer,m_lin_obj_factory);
1238  piro_params->sublist("NOX").sublist("Solver Options").set("User Defined Pre/Post Operator", ppo);
1239  }
1240 
1241  // override printing to use panzer ostream
1242  piro_params->sublist("NOX").sublist("Printing").set<Teuchos::RCP<std::ostream> >("Output Stream",global_data->os);
1243  piro_params->sublist("NOX").sublist("Printing").set<Teuchos::RCP<std::ostream> >("Error Stream",global_data->os);
1244  piro_params->sublist("NOX").sublist("Printing").set<int>("Output Processor",global_data->os->getOutputToRootOnly());
1245 
1246  // use the user specfied rythmos solver if they pass one in
1247  Teuchos::RCP<Piro::RythmosSolver<double> > piro_rythmos;
1248  if(rythmosSolver==Teuchos::null)
1249  piro_rythmos = Teuchos::rcp(new Piro::RythmosSolver<double>());
1250  else
1251  piro_rythmos = rythmosSolver;
1252 
1253  // if you are using explicit RK, make sure to wrap the ME in an explicit model evaluator decorator
1254  Teuchos::RCP<Thyra::ModelEvaluator<ScalarT> > rythmos_me = thyra_me;
1255  const std::string stepper_type = piro_params->sublist("Rythmos").get<std::string>("Stepper Type");
1256  if(stepper_type=="Explicit RK" || stepper_type=="Forward Euler") {
1257  const Teuchos::ParameterList & assembly_params = p.sublist("Assembly");
1258  bool lumpExplicitMass = assembly_params.get<bool>("Lump Explicit Mass");
1259  rythmos_me = Teuchos::rcp(new panzer::ExplicitModelEvaluator<ScalarT>(thyra_me,!useDynamicCoordinates_,lumpExplicitMass));
1260  }
1261 
1262  piro_rythmos->initialize(piro_params, rythmos_me, rythmos_observer_factory->buildRythmosObserver(m_mesh,m_global_indexer,m_lin_obj_factory));
1263 
1264  piro = piro_rythmos;
1265  }
1266  else {
1267  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
1268  "Error: Unknown Piro Solver : " << solver);
1269  }
1270  return piro;
1271  }
1272 
1273  template<typename ScalarT>
1274  Teuchos::RCP<panzer::ResponseLibrary<panzer::Traits> > ModelEvaluatorFactory<ScalarT>::getResponseLibrary()
1275  {
1276  TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::is_null(m_response_library), std::runtime_error,
1277  "Objects are not built yet! Please call buildObjects() member function.");
1278 
1279  return m_response_library;
1280  }
1281 
1282  template<typename ScalarT>
1283  const std::vector<Teuchos::RCP<panzer::PhysicsBlock> > & ModelEvaluatorFactory<ScalarT>::getPhysicsBlocks() const
1284  {
1285  TEUCHOS_TEST_FOR_EXCEPTION(m_physics_blocks.size()==0, std::runtime_error,
1286  "Objects are not built yet! Please call buildObjects() member function.");
1287 
1288  return m_physics_blocks;
1289  }
1290 
1291  template<typename ScalarT>
1292  Teuchos::RCP<panzer::FieldManagerBuilder>
1294  buildFieldManagerBuilder(const Teuchos::RCP<panzer::WorksetContainer> & wc,
1295  const std::vector<Teuchos::RCP<panzer::PhysicsBlock> >& physicsBlocks,
1296  const std::vector<panzer::BC> & bcs,
1297  const panzer::EquationSetFactory & eqset_factory,
1298  const panzer::BCStrategyFactory& bc_factory,
1301  const Teuchos::ParameterList& closure_models,
1302  const panzer::LinearObjFactory<panzer::Traits> & lo_factory,
1303  const Teuchos::ParameterList& user_data,
1304  bool writeGraph,const std::string & graphPrefix,
1305  bool write_field_managers,const std::string & field_manager_prefix) const
1306  {
1307  Teuchos::RCP<panzer::FieldManagerBuilder> fmb = Teuchos::rcp(new panzer::FieldManagerBuilder);
1308  fmb->setWorksetContainer(wc);
1309  fmb->setupVolumeFieldManagers(physicsBlocks,volume_cm_factory,closure_models,lo_factory,user_data);
1310  fmb->setupBCFieldManagers(bcs,physicsBlocks,eqset_factory,bc_cm_factory,bc_factory,closure_models,lo_factory,user_data);
1311 
1312  // Print Phalanx DAGs
1313  if (writeGraph){
1314  fmb->writeVolumeGraphvizDependencyFiles(graphPrefix, physicsBlocks);
1315  fmb->writeBCGraphvizDependencyFiles(graphPrefix);
1316  }
1317  if (write_field_managers){
1318  fmb->writeVolumeTextDependencyFiles(graphPrefix, physicsBlocks);
1319  fmb->writeBCTextDependencyFiles(field_manager_prefix);
1320  }
1321 
1322  return fmb;
1323  }
1324 
1325  template<typename ScalarT>
1326  Teuchos::RCP<Thyra::ModelEvaluator<double> >
1329  const Teuchos::RCP<Teuchos::ParameterList> & physics_block_plist,
1330  const Teuchos::RCP<const panzer::EquationSetFactory>& eqset_factory,
1331  const panzer::BCStrategyFactory & bc_factory,
1333  bool is_transient,bool is_explicit,
1334  const Teuchos::Ptr<const Teuchos::ParameterList> & bc_list,
1335  const Teuchos::RCP<Thyra::ModelEvaluator<ScalarT> > & physics_me_in) const
1336  {
1337  typedef panzer::ModelEvaluator<ScalarT> PanzerME;
1338 
1339  Teuchos::RCP<Thyra::ModelEvaluator<ScalarT> > physics_me = physics_me_in==Teuchos::null ? m_physics_me : physics_me_in;
1340 
1341  const Teuchos::ParameterList& p = *this->getParameterList();
1342 
1343  // build PhysicsBlocks
1344  std::vector<Teuchos::RCP<panzer::PhysicsBlock> > physicsBlocks;
1345  {
1346  const Teuchos::ParameterList & assembly_params = p.sublist("Assembly");
1347 
1348  // setup physical mappings and boundary conditions
1349  std::map<std::string,std::string> block_ids_to_physics_ids;
1350  panzer::buildBlockIdToPhysicsIdMap(block_ids_to_physics_ids, p.sublist("Block ID to Physics ID Mapping"));
1351 
1352  // build cell ( block id -> cell topology ) mapping
1353  std::map<std::string,Teuchos::RCP<const shards::CellTopology> > block_ids_to_cell_topo;
1354  for(std::map<std::string,std::string>::const_iterator itr=block_ids_to_physics_ids.begin();
1355  itr!=block_ids_to_physics_ids.end();++itr) {
1356  block_ids_to_cell_topo[itr->first] = m_mesh->getCellTopology(itr->first);
1357  TEUCHOS_ASSERT(block_ids_to_cell_topo[itr->first]!=Teuchos::null);
1358  }
1359 
1360  std::size_t workset_size = Teuchos::as<std::size_t>(assembly_params.get<int>("Workset Size"));
1361 
1362  panzer::buildPhysicsBlocks(block_ids_to_physics_ids,
1363  block_ids_to_cell_topo,
1364  physics_block_plist,
1365  assembly_params.get<int>("Default Integration Order"),
1366  workset_size,
1367  eqset_factory,
1368  m_global_data,
1369  is_transient,
1370  physicsBlocks);
1371  }
1372 
1373  // build FMB
1374  Teuchos::RCP<panzer::FieldManagerBuilder> fmb;
1375  {
1376  const Teuchos::ParameterList & user_data_params = p.sublist("User Data");
1377 
1378  bool write_dot_files = false;
1379  std::string prefix = "Cloned_";
1380 
1381  std::vector<panzer::BC> bcs;
1382  if(bc_list==Teuchos::null) {
1383  panzer::buildBCs(bcs, p.sublist("Boundary Conditions"), m_global_data);
1384  }
1385  else {
1386  panzer::buildBCs(bcs, *bc_list, m_global_data);
1387  }
1388 
1389  fmb = buildFieldManagerBuilder(// Teuchos::rcp_const_cast<panzer::WorksetContainer>(
1390  // m_response_library!=Teuchos::null ? m_response_library->getWorksetContainer()
1391  // : m_wkstContainer),
1392  m_wkstContainer,
1393  physicsBlocks,
1394  bcs,
1395  *eqset_factory,
1396  bc_factory,
1397  user_cm_factory,
1398  user_cm_factory,
1399  p.sublist("Closure Models"),
1400  *m_lin_obj_factory,
1401  user_data_params,
1402  write_dot_files,prefix,
1403  write_dot_files,prefix);
1404  }
1405 
1406  Teuchos::RCP<panzer::ResponseLibrary<panzer::Traits> > response_library
1407  = Teuchos::rcp(new panzer::ResponseLibrary<panzer::Traits>(m_wkstContainer,
1408  m_global_indexer,
1409  m_lin_obj_factory));
1410  // = Teuchos::rcp(new panzer::ResponseLibrary<panzer::Traits>(m_response_library->getWorksetContainer(),
1411  // m_response_library->getGlobalIndexer(),
1412  // m_response_library->getLinearObjFactory()));
1413 
1414  // using the FMB, build the model evaluator
1415  {
1416  // get nominal input values, make sure they match with internal me
1417  Thyra::ModelEvaluatorBase::InArgs<ScalarT> nomVals = physics_me->getNominalValues();
1418 
1419  // determine if this is a Epetra or Thyra ME
1420  Teuchos::RCP<Thyra::EpetraModelEvaluator> ep_thyra_me = Teuchos::rcp_dynamic_cast<Thyra::EpetraModelEvaluator>(physics_me);
1421  Teuchos::RCP<PanzerME> panzer_me = Teuchos::rcp_dynamic_cast<PanzerME>(physics_me);
1422  bool useThyra = true;
1423  if(ep_thyra_me!=Teuchos::null)
1424  useThyra = false;
1425 
1426  // get parameter names
1427  std::vector<Teuchos::RCP<Teuchos::Array<std::string> > > p_names(physics_me->Np());
1428  std::vector<Teuchos::RCP<Teuchos::Array<double> > > p_values(physics_me->Np());
1429  for(std::size_t i=0;i<p_names.size();i++) {
1430  p_names[i] = Teuchos::rcp(new Teuchos::Array<std::string>(*physics_me->get_p_names(i)));
1431  p_values[i] = Teuchos::rcp(new Teuchos::Array<double>(p_names[i]->size(),0.0));
1432  }
1433 
1434  Teuchos::RCP<Thyra::ModelEvaluatorDefaultBase<double> > thyra_me
1435  = buildPhysicsModelEvaluator(useThyra,
1436  fmb,
1437  response_library,
1438  m_lin_obj_factory,
1439  p_names,
1440  p_values,
1441  solverFactory,
1442  m_global_data,
1443  is_transient,
1444  nomVals.get_t());
1445 
1446  // set the nominal values...does this work???
1447  thyra_me->getNominalValues() = nomVals;
1448 
1449  // build an explicit model evaluator
1450  if(is_explicit) {
1451  const Teuchos::ParameterList & assembly_params = p.sublist("Assembly");
1452  bool lumpExplicitMass = assembly_params.get<bool>("Lump Explicit Mass");
1453  thyra_me = Teuchos::rcp(new panzer::ExplicitModelEvaluator<ScalarT>(thyra_me,!useDynamicCoordinates_,lumpExplicitMass));
1454  }
1455 
1456  return thyra_me;
1457  }
1458  }
1459 
1460  template<typename ScalarT>
1461  Teuchos::RCP<Thyra::ModelEvaluatorDefaultBase<double> >
1464  const Teuchos::RCP<panzer::FieldManagerBuilder> & fmb,
1465  const Teuchos::RCP<panzer::ResponseLibrary<panzer::Traits> > & rLibrary,
1466  const Teuchos::RCP<panzer::LinearObjFactory<panzer::Traits> > & lof,
1467  const std::vector<Teuchos::RCP<Teuchos::Array<std::string> > > & p_names,
1468  const std::vector<Teuchos::RCP<Teuchos::Array<double> > > & p_values,
1469  const Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<ScalarT> > & solverFactory,
1470  const Teuchos::RCP<panzer::GlobalData> & global_data,
1471  bool is_transient,double t_init) const
1472  {
1473  Teuchos::RCP<Thyra::ModelEvaluatorDefaultBase<double> > thyra_me;
1474  if(!buildThyraME) {
1475  Teuchos::RCP<panzer::ModelEvaluator_Epetra> ep_me
1476  = Teuchos::rcp(new panzer::ModelEvaluator_Epetra(fmb,rLibrary,lof, p_names,p_values, global_data, is_transient));
1477 
1478  if (is_transient)
1479  ep_me->set_t_init(t_init);
1480 
1481  // Build Thyra Model Evaluator
1482  thyra_me = Thyra::epetraModelEvaluator(ep_me,solverFactory);
1483  }
1484  else {
1485  thyra_me = Teuchos::rcp(new panzer::ModelEvaluator<double>
1486  (fmb,rLibrary,lof,p_names,p_values,solverFactory,global_data,is_transient,t_init));
1487  }
1488 
1489  return thyra_me;
1490  }
1491 
1492  template<typename ScalarT>
1494  getInitialTime(Teuchos::ParameterList& p,
1495  const panzer_stk::STK_Interface & mesh) const
1496  {
1497  Teuchos::ParameterList validPL;
1498  {
1499  Teuchos::setStringToIntegralParameter<int>(
1500  "Start Time Type",
1501  "From Input File",
1502  "Set the start time",
1503  Teuchos::tuple<std::string>("From Input File","From Exodus File"),
1504  &validPL
1505  );
1506 
1507  validPL.set<double>("Start Time",0.0);
1508  }
1509 
1510  p.validateParametersAndSetDefaults(validPL);
1511 
1512  std::string t_init_type = p.get<std::string>("Start Time Type");
1513  double t_init = 10.0;
1514 
1515  if (t_init_type == "From Input File")
1516  t_init = p.get<double>("Start Time");
1517 
1518  if (t_init_type == "From Exodus File")
1519  t_init = mesh.getInitialStateTime();
1520 
1521  return t_init;
1522  }
1523 
1524  // Setup STK response library for writing out the solution fields
1526  template<typename ScalarT>
1527  Teuchos::RCP<panzer::ResponseLibrary<panzer::Traits> > ModelEvaluatorFactory<ScalarT>::
1528  initializeSolnWriterResponseLibrary(const Teuchos::RCP<panzer::WorksetContainer> & wc,
1529  const Teuchos::RCP<const panzer::GlobalIndexer> & ugi,
1530  const Teuchos::RCP<const panzer::LinearObjFactory<panzer::Traits> > & lof,
1531  const Teuchos::RCP<panzer_stk::STK_Interface> & mesh) const
1532  {
1533  Teuchos::RCP<panzer::ResponseLibrary<panzer::Traits> > stkIOResponseLibrary
1534  = Teuchos::rcp(new panzer::ResponseLibrary<panzer::Traits>(wc,ugi,lof));
1535 
1536  std::vector<std::string> eBlocks;
1537  mesh->getElementBlockNames(eBlocks);
1538 
1540  builder.mesh = mesh;
1541  stkIOResponseLibrary->addResponse("Main Field Output",eBlocks,builder);
1542 
1543  return stkIOResponseLibrary;
1544  }
1545 
1546  template<typename ScalarT>
1549  const std::vector<Teuchos::RCP<panzer::PhysicsBlock> > & physicsBlocks,
1551  const Teuchos::ParameterList & closure_models,
1552  int workset_size, Teuchos::ParameterList & user_data) const
1553  {
1554  user_data.set<int>("Workset Size",workset_size);
1555  rl.buildResponseEvaluators(physicsBlocks, cm_factory, closure_models, user_data);
1556  }
1557 
1558  template<typename ScalarT>
1559  Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<double> > ModelEvaluatorFactory<ScalarT>::
1560  buildLOWSFactory(bool blockedAssembly,
1561  const Teuchos::RCP<const panzer::GlobalIndexer> & globalIndexer,
1562  const Teuchos::RCP<panzer::ConnManager> & conn_manager,
1563  const Teuchos::RCP<panzer_stk::STK_Interface> & mesh,
1564  const Teuchos::RCP<const Teuchos::MpiComm<int> > & mpi_comm
1565  #ifdef PANZER_HAVE_TEKO
1566  , const Teuchos::RCP<Teko::RequestHandler> & reqHandler
1567  #endif
1568  ) const
1569  {
1570  const Teuchos::ParameterList & p = *this->getParameterList();
1571  const Teuchos::ParameterList & solncntl_params = p.sublist("Solution Control");
1572 
1573  // Build stratimikos solver (note that this is a hard coded path to linear solver options in nox list!)
1574  Teuchos::RCP<Teuchos::ParameterList> strat_params
1575  = Teuchos::rcp(new Teuchos::ParameterList(solncntl_params.sublist("NOX").sublist("Direction").
1576  sublist("Newton").sublist("Stratimikos Linear Solver").sublist("Stratimikos")));
1577 
1578  bool writeCoordinates = false;
1579  if(p.sublist("Options").isType<bool>("Write Coordinates"))
1580  writeCoordinates = p.sublist("Options").get<bool>("Write Coordinates");
1581 
1582  bool writeTopo = false;
1583  if(p.sublist("Options").isType<bool>("Write Topology"))
1584  writeTopo = p.sublist("Options").get<bool>("Write Topology");
1585 
1586 
1588  blockedAssembly,globalIndexer,conn_manager,
1589  Teuchos::as<int>(mesh->getDimension()), mpi_comm, strat_params,
1590  #ifdef PANZER_HAVE_TEKO
1591  reqHandler,
1592  #endif
1593  writeCoordinates,
1594  writeTopo
1595  );
1596  }
1597 
1598  template<typename ScalarT>
1601  const bool write_graphviz_file,
1602  const std::string& graphviz_file_prefix)
1603  {
1604  typedef panzer::ModelEvaluator<double> PanzerME;
1605 
1606  Teuchos::ParameterList & p = *this->getNonconstParameterList();
1607  Teuchos::ParameterList & user_data = p.sublist("User Data");
1608  Teuchos::ParameterList & closure_models = p.sublist("Closure Models");
1609 
1610  // uninitialize the thyra model evaluator, its respone counts are wrong!
1611  Teuchos::RCP<Thyra::EpetraModelEvaluator> thyra_me = Teuchos::rcp_dynamic_cast<Thyra::EpetraModelEvaluator>(m_physics_me);
1612  Teuchos::RCP<PanzerME> panzer_me = Teuchos::rcp_dynamic_cast<PanzerME>(m_physics_me);
1613 
1614  if(thyra_me!=Teuchos::null && panzer_me==Teuchos::null) {
1615  Teuchos::RCP<const EpetraExt::ModelEvaluator> const_ep_me;
1616  Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<double> > solveFactory;
1617  thyra_me->uninitialize(&const_ep_me,&solveFactory); // this seems dangerous!
1618 
1619  // I don't need no const-ness!
1620  Teuchos::RCP<EpetraExt::ModelEvaluator> ep_me = Teuchos::rcp_const_cast<EpetraExt::ModelEvaluator>(const_ep_me);
1621  Teuchos::RCP<panzer::ModelEvaluator_Epetra> ep_panzer_me = Teuchos::rcp_dynamic_cast<panzer::ModelEvaluator_Epetra>(ep_me);
1622 
1623  ep_panzer_me->buildResponses(m_physics_blocks,*m_eqset_factory,cm_factory,closure_models,user_data,write_graphviz_file,graphviz_file_prefix);
1624 
1625  // reinitialize the thyra model evaluator, now with the correct responses
1626  thyra_me->initialize(ep_me,solveFactory);
1627 
1628  return;
1629  }
1630  else if(panzer_me!=Teuchos::null && thyra_me==Teuchos::null) {
1631  panzer_me->buildResponses(m_physics_blocks,*m_eqset_factory,cm_factory,closure_models,user_data,write_graphviz_file,graphviz_file_prefix);
1632 
1633  return;
1634  }
1635 
1636  TEUCHOS_ASSERT(false);
1637  }
1638 }
1639 
1640 #endif
void TokensToInts(std::vector< int > &values, const std::vector< std::string > &tokens)
Turn a vector of tokens into a vector of ints.
void setRythmosObserverFactory(const Teuchos::RCP< const panzer_stk::RythmosObserverFactory > &rythmos_observer_factory)
Interface for constructing a BCStrategy_TemplateManager.
virtual void completeMeshConstruction(STK_Interface &mesh, stk::ParallelMachine parallelMach) const =0
Allocates and initializes an equation set template manager.
void finalizeMeshConstruction(const STK_MeshFactory &mesh_factory, const std::vector< Teuchos::RCP< panzer::PhysicsBlock > > &physicsBlocks, const Teuchos::MpiComm< int > mpi_comm, STK_Interface &mesh) const
Teuchos::RCP< Thyra::ModelEvaluator< ScalarT > > getPhysicsModelEvaluator()
Teuchos::RCP< Thyra::ModelEvaluator< double > > cloneWithNewPhysicsBlocks(const Teuchos::RCP< Thyra::LinearOpWithSolveFactoryBase< ScalarT > > &solverFactory, const Teuchos::RCP< Teuchos::ParameterList > &physics_block_plist, const Teuchos::RCP< const panzer::EquationSetFactory > &eqset_factory, const panzer::BCStrategyFactory &bc_factory, const panzer::ClosureModelFactory_TemplateManager< panzer::Traits > &user_cm_factory, bool is_transient, bool is_explicit, const Teuchos::Ptr< const Teuchos::ParameterList > &bc_list=Teuchos::null, const Teuchos::RCP< Thyra::ModelEvaluator< ScalarT > > &physics_me=Teuchos::null) const
void writeInitialConditions(const Thyra::ModelEvaluator< ScalarT > &model, const std::vector< Teuchos::RCP< panzer::PhysicsBlock > > &physicsBlocks, const Teuchos::RCP< panzer::WorksetContainer > &wc, const Teuchos::RCP< const panzer::GlobalIndexer > &ugi, const Teuchos::RCP< const panzer::LinearObjFactory< panzer::Traits > > &lof, const Teuchos::RCP< panzer_stk::STK_Interface > &mesh, const panzer::ClosureModelFactory_TemplateManager< panzer::Traits > &cm_factory, const Teuchos::ParameterList &closure_model_pl, const Teuchos::ParameterList &user_data_pl, int workset_size) const
Write the initial conditions to exodus. Note that this is entirely self contained.
static bool requiresBlocking(const std::string &fieldorder)
virtual Teuchos::RCP< panzer::GlobalIndexer > buildGlobalIndexer(const Teuchos::RCP< const Teuchos::OpaqueWrapper< MPI_Comm > > &mpiComm, const std::vector< Teuchos::RCP< panzer::PhysicsBlock > > &physicsBlocks, const Teuchos::RCP< ConnManager > &connMngr, const std::string &fieldOrder="") const
void buildBlockIdToPhysicsIdMap(std::map< std::string, std::string > &b_to_p, const Teuchos::ParameterList &p)
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
virtual void set_x_th(const Teuchos::RCP< Thyra::VectorBase< ScalarT > > &in)=0
void addSolutionField(const std::string &fieldName, const std::string &blockId)
void setParameterList(Teuchos::RCP< Teuchos::ParameterList > const &paramList)
Teuchos::RCP< panzer::LinearObjContainer > ghostedContainer_
Class that provides access to worksets on each element block and side set.
Teuchos::RCP< panzer::LinearObjContainer > container_
Teuchos::RCP< panzer::FieldManagerBuilder > buildFieldManagerBuilder(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 panzer::LinearObjFactory< panzer::Traits > &lo_factory, const Teuchos::ParameterList &user_data, bool writeGraph, const std::string &graphPrefix, bool write_field_managers, const std::string &field_manager_prefix) const
unsigned getDimension() const
get the dimension
void buildResponseEvaluators(const std::vector< Teuchos::RCP< panzer::PhysicsBlock > > &physicsBlocks, 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< Thyra::ModelEvaluator< ScalarT > > buildResponseOnlyModelEvaluator(const Teuchos::RCP< Thyra::ModelEvaluator< ScalarT > > &thyra_me, const Teuchos::RCP< panzer::GlobalData > &global_data, const Teuchos::RCP< Piro::RythmosSolver< ScalarT > > rythmosSolver=Teuchos::null, const Teuchos::Ptr< const panzer_stk::NOXObserverFactory > &in_nox_observer_factory=Teuchos::null, const Teuchos::Ptr< const panzer_stk::RythmosObserverFactory > &in_rythmos_observer_factory=Teuchos::null)
void buildResponses(const std::vector< Teuchos::RCP< panzer::PhysicsBlock > > &physicsBlocks, const panzer::EquationSetFactory &eqset_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="")
void addUserFieldsToMesh(panzer_stk::STK_Interface &mesh, const Teuchos::ParameterList &output_list) const
Add the user fields specified by output_list to the mesh.
Teuchos::RCP< Thyra::ModelEvaluator< ScalarT > > getResponseOnlyModelEvaluator()
Teuchos::RCP< Thyra::ModelEvaluatorDefaultBase< double > > buildPhysicsModelEvaluator(bool buildThyraME, 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< Thyra::LinearOpWithSolveFactoryBase< ScalarT > > &solverFactory, const Teuchos::RCP< panzer::GlobalData > &global_data, bool is_transient, double t_init) const
void setNOXObserverFactory(const Teuchos::RCP< const panzer_stk::NOXObserverFactory > &nox_observer_factory)
Teuchos::RCP< STK_MeshFactory > buildSTKMeshFactory(const Teuchos::ParameterList &mesh_params) const
build STK mesh factory from a mesh parameter list
Teuchos::RCP< Thyra::LinearOpWithSolveFactoryBase< double > > buildLOWSFactory(bool blockedAssembly, const Teuchos::RCP< const panzer::GlobalIndexer > &globalIndexer, const Teuchos::RCP< panzer::ConnManager > &conn_manager, const Teuchos::RCP< panzer_stk::STK_Interface > &mesh, const Teuchos::RCP< const Teuchos::MpiComm< int > > &mpi_comm) const
void buildObjects(const Teuchos::RCP< const Teuchos::Comm< int > > &comm, const Teuchos::RCP< panzer::GlobalData > &global_data, const Teuchos::RCP< const panzer::EquationSetFactory > &eqset_factory, const panzer::BCStrategyFactory &bc_factory, const panzer::ClosureModelFactory_TemplateManager< panzer::Traits > &cm_factory, bool meConstructionOn=true)
Builds the model evaluators for a panzer assembly.
virtual Teuchos::RCP< LinearObjContainer > buildLinearObjContainer() const =0
void setBlockWeight(const std::string &blockId, double weight)
void setupInitialConditions(Thyra::ModelEvaluator< ScalarT > &model, panzer::WorksetContainer &wkstContainer, const std::vector< Teuchos::RCP< panzer::PhysicsBlock > > &physicsBlocks, const panzer::ClosureModelFactory_TemplateManager< panzer::Traits > &cm_factory, const panzer::LinearObjFactory< panzer::Traits > &lof, const Teuchos::ParameterList &closure_pl, const Teuchos::ParameterList &initial_cond_pl, const Teuchos::ParameterList &user_data_pl, bool write_dot_files, const std::string &dot_file_prefix) const
Setup the initial conditions in a model evaluator. Note that this is entirely self contained...
Teuchos::RCP< Thyra::LinearOpWithSolveFactoryBase< double > > buildLOWSFactory(bool blockedAssembly, const Teuchos::RCP< const panzer::GlobalIndexer > &globalIndexer, const Teuchos::RCP< panzer_stk::STKConnManager > &stkConn_manager, int spatialDim, const Teuchos::RCP< const Teuchos::MpiComm< int > > &mpi_comm, const Teuchos::RCP< Teuchos::ParameterList > &strat_params, bool writeCoordinates, bool writeTopo, const Teuchos::RCP< const panzer::GlobalIndexer > &auxGlobalIndexer, bool useCoordinates)
double getInitialTime(Teuchos::ParameterList &transient_ic_params, const panzer_stk::STK_Interface &mesh) const
Gets the initial time from either the input parameter list or an exodus file.
virtual void readVector(const std::string &identifier, LinearObjContainer &loc, int id) const =0
void addMeshCoordFields(const std::string &blockId, const std::vector< std::string > &coordField, const std::string &dispPrefix)
virtual Teuchos::RCP< panzer::GlobalIndexer > buildGlobalIndexer(const Teuchos::RCP< const Teuchos::OpaqueWrapper< MPI_Comm > > &mpiComm, const std::vector< Teuchos::RCP< panzer::PhysicsBlock > > &physicsBlocks, const Teuchos::RCP< ConnManager > &connMngr, const std::string &fieldOrder="") const
void setupInitialConditionFieldManagers(WorksetContainer &wkstContainer, const std::vector< Teuchos::RCP< panzer::PhysicsBlock > > &physicsBlocks, const panzer::ClosureModelFactory_TemplateManager< panzer::Traits > &cm_factory, const Teuchos::ParameterList &ic_block_closure_models, const panzer::LinearObjFactory< panzer::Traits > &lo_factory, const Teuchos::ParameterList &user_data, const bool write_graphviz_file, const std::string &graphviz_file_prefix, std::map< std::string, Teuchos::RCP< PHX::FieldManager< panzer::Traits > > > &phx_ic_field_managers)
Builds PHX::FieldManager objects for inital conditions and registers evaluators.
void buildResponses(const panzer::ClosureModelFactory_TemplateManager< panzer::Traits > &cm_factory, const bool write_graphviz_file=false, const std::string &graphviz_file_prefix="")
const std::vector< Teuchos::RCP< panzer::PhysicsBlock > > & getPhysicsBlocks() const
void finalizeSolnWriterResponseLibrary(panzer::ResponseLibrary< panzer::Traits > &rl, const std::vector< Teuchos::RCP< panzer::PhysicsBlock > > &physicsBlocks, const panzer::ClosureModelFactory_TemplateManager< panzer::Traits > &cm_factory, const Teuchos::ParameterList &closure_models, int workset_size, Teuchos::ParameterList &user_data) const
std::pair< std::string, Teuchos::RCP< panzer::PureBasis > > StrPureBasisPair
Teuchos::RCP< panzer::ResponseLibrary< panzer::Traits > > initializeSolnWriterResponseLibrary(const Teuchos::RCP< panzer::WorksetContainer > &wc, const Teuchos::RCP< const panzer::GlobalIndexer > &ugi, const Teuchos::RCP< const panzer::LinearObjFactory< panzer::Traits > > &lof, const Teuchos::RCP< panzer_stk::STK_Interface > &mesh) const
WorksetDescriptor blockDescriptor(const std::string &eBlock)
void StringTokenizer(std::vector< std::string > &tokens, const std::string &str, const std::string delimiters, bool trim)
Tokenize a string, put tokens in a vector.
void addCellField(const std::string &fieldName, const std::string &blockId)
void evaluateInitialCondition(WorksetContainer &wkstContainer, const std::map< std::string, Teuchos::RCP< PHX::FieldManager< panzer::Traits > > > &phx_ic_field_managers, Teuchos::RCP< panzer::LinearObjContainer > loc, const panzer::LinearObjFactory< panzer::Traits > &lo_factory, const double time_stamp)
void registerScalarParameter(const std::string name, panzer::ParamLib &pl, double realValue)
void setUserWorksetFactory(Teuchos::RCP< panzer_stk::WorksetFactory > &user_wkst_factory)
Set user defined workset factory.
Teuchos::RCP< panzer::ResponseLibrary< panzer::Traits > > getResponseLibrary()
std::string printUGILoadBalancingInformation(const GlobalIndexer &ugi)