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