Panzer  Version of the Day
Panzer_STK_SetupLOWSFactory_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_SetupLOWSFactory_impl_hpp__
44 #define __Panzer_STK_SetupLOWSFactory_impl_hpp__
45 
46 #include "PanzerAdaptersSTK_config.hpp"
49 
50 #include "Teuchos_AbstractFactoryStd.hpp"
51 
52 #include "Stratimikos_DefaultLinearSolverBuilder.hpp"
53 
54 #include "Epetra_MpiComm.h"
55 #include "Epetra_Vector.h"
56 #include "EpetraExt_VectorOut.h"
57 
58 #include "Tpetra_Map.hpp"
59 #include "Tpetra_MultiVector.hpp"
60 
61 #ifdef PANZER_HAVE_TEKO
62 #include "Teko_StratimikosFactory.hpp"
63 #endif
64 
65 #ifdef PANZER_HAVE_MUELU
66 #include <Thyra_MueLuPreconditionerFactory.hpp>
67 #include "Stratimikos_MueLuHelpers.hpp"
68 #include "MatrixMarket_Tpetra.hpp"
69 #endif
70 
71 #ifdef PANZER_HAVE_IFPACK2
72 #include <Thyra_Ifpack2PreconditionerFactory.hpp>
73 #endif
74 
75 namespace panzer_stk {
76 
77 namespace {
78 
79  bool
80  determineCoordinateField(const panzer::UniqueGlobalIndexerBase & globalIndexer,std::string & fieldName)
81  {
82  std::vector<std::string> elementBlocks;
83  globalIndexer.getElementBlockIds(elementBlocks);
84 
85  // grab fields for first block
86  std::set<int> runningFields;
87  {
88  const std::vector<int> & fields = globalIndexer.getBlockFieldNumbers(elementBlocks[0]);
89  runningFields.insert(fields.begin(),fields.end());
90  }
91 
92  // grab fields for first block
93  for(std::size_t i=1;i<elementBlocks.size();i++) {
94  const std::vector<int> & fields = globalIndexer.getBlockFieldNumbers(elementBlocks[i]);
95 
96  std::set<int> currentFields(runningFields);
97  runningFields.clear();
98  std::set_intersection(fields.begin(),fields.end(),
99  currentFields.begin(),currentFields.end(),
100  std::inserter(runningFields,runningFields.begin()));
101  }
102 
103  if(runningFields.size()<1)
104  return false;
105 
106  fieldName = globalIndexer.getFieldString(*runningFields.begin());
107  return true;
108  }
109 
110 #ifdef PANZER_HAVE_FEI
111  template<typename GO>
112  void
113  fillFieldPatternMap(const panzer::DOFManagerFEI<int,GO> & globalIndexer,
114  const std::string & fieldName,
115  std::map<std::string,Teuchos::RCP<const panzer::Intrepid2FieldPattern> > & fieldPatterns)
116  {
117  std::vector<std::string> elementBlocks;
118  globalIndexer.getElementBlockIds(elementBlocks);
119 
120  for(std::size_t e=0;e<elementBlocks.size();e++) {
121  std::string blockId = elementBlocks[e];
122 
123  if(globalIndexer.fieldInBlock(fieldName,blockId))
124  fieldPatterns[blockId] =
125  Teuchos::rcp_dynamic_cast<const panzer::Intrepid2FieldPattern>(globalIndexer.getFieldPattern(blockId,fieldName),true);
126  }
127  }
128 #endif
129 
130  template<typename GO>
131  void
132  fillFieldPatternMap(const panzer::DOFManager<int,GO> & globalIndexer,
133  const std::string & fieldName,
134  std::map<std::string,Teuchos::RCP<const panzer::Intrepid2FieldPattern> > & fieldPatterns)
135  {
136  std::vector<std::string> elementBlocks;
137  globalIndexer.getElementBlockIds(elementBlocks);
138 
139  for(std::size_t e=0;e<elementBlocks.size();e++) {
140  std::string blockId = elementBlocks[e];
141 
142  if(globalIndexer.fieldInBlock(fieldName,blockId))
143  fieldPatterns[blockId] =
144  Teuchos::rcp_dynamic_cast<const panzer::Intrepid2FieldPattern>(globalIndexer.getFieldPattern(blockId,fieldName),true);
145  }
146  }
147 
148  void
149  fillFieldPatternMap(const panzer::UniqueGlobalIndexerBase & globalIndexer,
150  const std::string & fieldName,
151  std::map<std::string,Teuchos::RCP<const panzer::Intrepid2FieldPattern> > & fieldPatterns)
152  {
153  using Teuchos::Ptr;
154  using Teuchos::ptrFromRef;
155  using Teuchos::ptr_dynamic_cast;
156  using panzer::DOFManager;
157 #ifdef PANZER_HAVE_FEI
158  using panzer::DOFManagerFEI;
159 #endif
160 
161  // first standard dof manager
162  {
163  Ptr<const DOFManager<int,int> > dofManager = ptr_dynamic_cast<const DOFManager<int,int> >(ptrFromRef(globalIndexer));
164 
165  if(dofManager!=Teuchos::null) {
166  fillFieldPatternMap(*dofManager,fieldName,fieldPatterns);
167  return;
168  }
169  }
170  {
171  Ptr<const DOFManager<int,panzer::Ordinal64> > dofManager = ptr_dynamic_cast<const DOFManager<int,panzer::Ordinal64> >(ptrFromRef(globalIndexer));
172 
173  if(dofManager!=Teuchos::null) {
174  fillFieldPatternMap(*dofManager,fieldName,fieldPatterns);
175  return;
176  }
177  }
178 
179 #ifdef PANZER_HAVE_FEI
180  // now FEI dof manager
181  {
182  Ptr<const DOFManagerFEI<int,int> > dofManager = ptr_dynamic_cast<const DOFManagerFEI<int,int> >(ptrFromRef(globalIndexer));
183 
184  if(dofManager!=Teuchos::null) {
185  fillFieldPatternMap(*dofManager,fieldName,fieldPatterns);
186  return;
187  }
188  }
189  {
190  Ptr<const DOFManagerFEI<int,panzer::Ordinal64> > dofManager = ptr_dynamic_cast<const DOFManagerFEI<int,panzer::Ordinal64> >(ptrFromRef(globalIndexer));
191 
192  if(dofManager!=Teuchos::null) {
193  fillFieldPatternMap(*dofManager,fieldName,fieldPatterns);
194  return;
195  }
196  }
197 #endif
198  }
199 }
200 
201  template<typename GO>
203  buildLOWSFactory(bool blockedAssembly,
205  const Teuchos::RCP<panzer_stk::STKConnManager<GO> > & stkConn_manager,
206  int spatialDim,
207  const Teuchos::RCP<const Teuchos::MpiComm<int> > & mpi_comm,
208  const Teuchos::RCP<Teuchos::ParameterList> & strat_params,
209  #ifdef PANZER_HAVE_TEKO
210  const Teuchos::RCP<Teko::RequestHandler> & reqHandler,
211  #endif
212  bool writeCoordinates,
213  bool writeTopo
214  )
215  {
216  Stratimikos::DefaultLinearSolverBuilder linearSolverBuilder;
217 
218  // Note if you want to use new solvers within Teko they have to be added to the solver builer
219  // before teko is added. This is because Teko steals its defaults from the solver its being injected
220  // into!
221 
222  #ifdef PANZER_HAVE_MUELU
223  {
224  // TAW: the following is probably not optimal but it corresponds to what have been there before...
225  Stratimikos::enableMueLu(linearSolverBuilder,"MueLu");
226  Stratimikos::enableMueLu<int,panzer::Ordinal64,panzer::TpetraNodeType>(linearSolverBuilder,"MueLu-Tpetra");
227  }
228  #endif // MUELU
229  #ifdef PANZER_HAVE_IFPACK2
230  {
231  typedef Thyra::PreconditionerFactoryBase<double> Base;
232  typedef Thyra::Ifpack2PreconditionerFactory<Tpetra::CrsMatrix<double, int, panzer::Ordinal64,panzer::TpetraNodeType> > Impl;
233 
234  linearSolverBuilder.setPreconditioningStrategyFactory(Teuchos::abstractFactoryStd<Base, Impl>(), "Ifpack2");
235  }
236  #endif // MUELU
237 
238 
239  #ifdef PANZER_HAVE_TEKO
240  Teuchos::RCP<Teko::RequestHandler> reqHandler_local = reqHandler;
241 
242  if(!blockedAssembly) {
243 
244  std::string fieldName;
245 
246  // try to set request handler from member variable
247  if(reqHandler_local==Teuchos::null)
248  reqHandler_local = Teuchos::rcp(new Teko::RequestHandler);
249 
250  // add in the coordinate parameter list callback handler
251  if(determineCoordinateField(*globalIndexer,fieldName)) {
252  std::map<std::string,Teuchos::RCP<const panzer::Intrepid2FieldPattern> > fieldPatterns;
253  fillFieldPatternMap(*globalIndexer,fieldName,fieldPatterns);
254 
256  panzer_stk::ParameterListCallback<int,GO>(fieldName,fieldPatterns,stkConn_manager,
257  Teuchos::rcp_dynamic_cast<const panzer::UniqueGlobalIndexer<int,GO> >(globalIndexer)));
258  reqHandler_local->addRequestCallback(callback);
259 
260  if(writeCoordinates) {
261  // force parameterlistcallback to build coordinates
262  callback->preRequest(Teko::RequestMesg(Teuchos::rcp(new Teuchos::ParameterList())));
263 
264  // extract coordinate vectors
265  const std::vector<double> & xcoords = callback->getXCoordsVector();
266  const std::vector<double> & ycoords = callback->getYCoordsVector();
267  const std::vector<double> & zcoords = callback->getZCoordsVector();
268 
269  // use epetra to write coordinates to matrix market files
270  Epetra_MpiComm ep_comm(*mpi_comm->getRawMpiComm()); // this is OK access to RawMpiComm becase its declared on the stack?
271  // and all users of this object are on the stack (within scope of mpi_comm
272  Epetra_Map map(-1,xcoords.size(),0,ep_comm);
273 
275  switch(spatialDim) {
276  case 3:
277  vec = Teuchos::rcp(new Epetra_Vector(Copy,map,const_cast<double *>(&zcoords[0])));
278  EpetraExt::VectorToMatrixMarketFile("zcoords.mm",*vec);
279  case 2:
280  vec = Teuchos::rcp(new Epetra_Vector(Copy,map,const_cast<double *>(&ycoords[0])));
281  EpetraExt::VectorToMatrixMarketFile("ycoords.mm",*vec);
282  case 1:
283  vec = Teuchos::rcp(new Epetra_Vector(Copy,map,const_cast<double *>(&xcoords[0])));
284  EpetraExt::VectorToMatrixMarketFile("xcoords.mm",*vec);
285  break;
286  default:
287  TEUCHOS_ASSERT(false);
288  }
289  }
290 
291  #ifdef PANZER_HAVE_MUELU
292  if(Teuchos::rcp_dynamic_cast<const panzer::UniqueGlobalIndexer<int,panzer::Ordinal64> >(globalIndexer)!=Teuchos::null) {
293  if(!writeCoordinates)
294  callback->preRequest(Teko::RequestMesg(Teuchos::rcp(new Teuchos::ParameterList())));
295 
296  typedef Tpetra::Map<int,panzer::Ordinal64,panzer::TpetraNodeType> Map;
297  typedef Tpetra::MultiVector<double,int,panzer::Ordinal64,panzer::TpetraNodeType> MV;
298 
299  // extract coordinate vectors and modify strat_params to include coordinate vectors
300  unsigned dim = Teuchos::as<unsigned>(spatialDim);
301  Teuchos::RCP<MV> coords;
302  for(unsigned d=0;d<dim;d++) {
303  const std::vector<double> & coord = callback->getCoordsVector(d);
304 
305  // no coords vector has been build yet, build one
306  if(coords==Teuchos::null) {
307  if(globalIndexer->getNumFields()==1) {
309  = Teuchos::rcp_dynamic_cast<const panzer::UniqueGlobalIndexer<int,panzer::Ordinal64> >(globalIndexer);
310  std::vector<panzer::Ordinal64> ownedIndices;
311  ugi->getOwnedIndices(ownedIndices);
312  Teuchos::RCP<const Map> coords_map = Teuchos::rcp(new Map(Teuchos::OrdinalTraits<panzer::Ordinal64>::invalid(),ownedIndices,0,mpi_comm));
313  coords = Teuchos::rcp(new MV(coords_map,dim));
314  }
315  else {
316  Teuchos::RCP<const Map> coords_map = Teuchos::rcp(new Map(Teuchos::OrdinalTraits<panzer::Ordinal64>::invalid(),coord.size(),0,mpi_comm));
317  coords = Teuchos::rcp(new MV(coords_map,dim));
318  }
319  }
320 
321  // sanity check the size
322  TEUCHOS_ASSERT(coords->getLocalLength()==coord.size());
323 
324  // fill appropriate coords vector
325  Teuchos::ArrayRCP<double> dest = coords->getDataNonConst(d);
326  for(std::size_t i=0;i<coord.size();i++)
327  dest[i] = coord[i];
328  }
329 
330  // inject coordinates into parameter list
331  Teuchos::ParameterList & muelu_params = strat_params->sublist("Preconditioner Types").sublist("MueLu-Tpetra");
332  muelu_params.set<Teuchos::RCP<MV> >("Coordinates",coords);
333  }
334  #endif
335  }
336  // else write_out_the_mesg("Warning: No unique field determines the coordinates, coordinates unavailable!")
337 
338  Teko::addTekoToStratimikosBuilder(linearSolverBuilder,reqHandler_local);
339  }
340  else {
341  // try to set request handler from member variable
342  if(reqHandler_local==Teuchos::null)
343  reqHandler_local = Teuchos::rcp(new Teko::RequestHandler);
344 
345  std::string fieldName;
346  if(determineCoordinateField(*globalIndexer,fieldName)) {
348  Teuchos::rcp_dynamic_cast<const panzer::BlockedDOFManager<int,GO> >(globalIndexer);
350  Teuchos::rcp(new panzer_stk::ParameterListCallbackBlocked<int,GO>(stkConn_manager,blkDofs));
351  reqHandler_local->addRequestCallback(callback);
352  }
353 
354  Teko::addTekoToStratimikosBuilder(linearSolverBuilder,reqHandler_local);
355 
356  if(writeCoordinates) {
358  Teuchos::rcp_dynamic_cast<const panzer::BlockedDOFManager<int,GO> >(globalIndexer);
359 
360  // loop over blocks
361  const std::vector<Teuchos::RCP<panzer::UniqueGlobalIndexer<int,GO> > > & dofVec
362  = blkDofs->getFieldDOFManagers();
363  for(std::size_t i=0;i<dofVec.size();i++) {
364  std::string fieldName;
365 
366  // add in the coordinate parameter list callback handler
367  TEUCHOS_ASSERT(determineCoordinateField(*dofVec[i],fieldName));
368 
369  std::map<std::string,Teuchos::RCP<const panzer::Intrepid2FieldPattern> > fieldPatterns;
370  fillFieldPatternMap(*dofVec[i],fieldName,fieldPatterns);
371  panzer_stk::ParameterListCallback<int,GO> plCall(fieldName,fieldPatterns,stkConn_manager,dofVec[i]);
372  plCall.buildArrayToVector();
373  plCall.buildCoordinates();
374 
375  // extract coordinate vectors
376  const std::vector<double> & xcoords = plCall.getXCoordsVector();
377  const std::vector<double> & ycoords = plCall.getYCoordsVector();
378  const std::vector<double> & zcoords = plCall.getZCoordsVector();
379 
380  // use epetra to write coordinates to matrix market files
381  Epetra_MpiComm ep_comm(*mpi_comm->getRawMpiComm()); // this is OK access to RawMpiComm becase its declared on the stack?
382  // and all users of this object are on the stack (within scope of mpi_comm
383  Epetra_Map map(-1,xcoords.size(),0,ep_comm);
384 
386  switch(spatialDim) {
387  case 3:
388  vec = Teuchos::rcp(new Epetra_Vector(Copy,map,const_cast<double *>(&zcoords[0])));
389  EpetraExt::VectorToMatrixMarketFile((fieldName+"_zcoords.mm").c_str(),*vec);
390  case 2:
391  vec = Teuchos::rcp(new Epetra_Vector(Copy,map,const_cast<double *>(&ycoords[0])));
392  EpetraExt::VectorToMatrixMarketFile((fieldName+"_ycoords.mm").c_str(),*vec);
393  case 1:
394  vec = Teuchos::rcp(new Epetra_Vector(Copy,map,const_cast<double *>(&xcoords[0])));
395  EpetraExt::VectorToMatrixMarketFile((fieldName+"_xcoords.mm").c_str(),*vec);
396  break;
397  default:
398  TEUCHOS_ASSERT(false);
399  }
400  }
401  }
402 
403  if(writeTopo) {
405  Teuchos::rcp_dynamic_cast<const panzer::BlockedDOFManager<int,GO> >(globalIndexer);
406 
407  writeTopology(*blkDofs);
408  }
409  }
410  #endif
411 
412  linearSolverBuilder.setParameterList(strat_params);
413  Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<double> > lowsFactory = createLinearSolveStrategy(linearSolverBuilder);
414 
415  return lowsFactory;
416  }
417 
418  template<typename GO>
419  void
421  {
422  using Teuchos::RCP;
423 
424  // loop over each field block
425  const std::vector<RCP<panzer::UniqueGlobalIndexer<int,GO> > > & blk_dofMngrs = blkDofs.getFieldDOFManagers();
426  for(std::size_t b=0;b<blk_dofMngrs.size();b++) {
427 #ifdef PANZER_HAVE_FEI
428  RCP<panzer::DOFManagerFEI<int,GO> > dofMngr = Teuchos::rcp_dynamic_cast<panzer::DOFManagerFEI<int,GO> >(blk_dofMngrs[b],true);
429 
430  std::vector<std::string> eBlocks;
431  dofMngr->getElementBlockIds(eBlocks);
432 
433  // build file name
434  std::stringstream fileName;
435  fileName << "elements_" << b;
436  std::ofstream file(fileName.str().c_str());
437 
438  // loop over each element block, write out topology
439  for(std::size_t e=0;e<eBlocks.size();e++)
440  writeTopology(*dofMngr,eBlocks[e],file);
441 #else
442  TEUCHOS_ASSERT(false);
443 #endif
444  }
445  }
446 
447 #ifdef PANZER_HAVE_FEI
448  template <typename GO>
449  void
450  writeTopology(const panzer::DOFManagerFEI<int,GO> & dofs,const std::string & block,std::ostream & os)
451  {
452  std::vector<std::string> fields(dofs.getElementBlockGIDCount(block));
453 
454  const std::set<int> & fieldIds = dofs.getFields(block);
455  for(std::set<int>::const_iterator itr=fieldIds.begin();itr!=fieldIds.end();++itr) {
456  std::string field = dofs.getFieldString(*itr);
457 
458  // get the layout of each field
459  const std::vector<int> & fieldOffsets = dofs.getGIDFieldOffsets(block,*itr);
460  for(std::size_t f=0;f<fieldOffsets.size();f++)
461  fields[fieldOffsets[f]] = field;
462 
463  }
464 
465  // print the layout of the full pattern
466  os << "#" << std::endl;
467  os << "# Element Block \"" << block << "\"" << std::endl;
468  os << "# field pattern = [ " << fields[0];
469  for(std::size_t f=1;f<fields.size();f++)
470  os << ", " << fields[f];
471  os << " ]" << std::endl;
472  os << "#" << std::endl;
473 
474  const std::vector<int> & elements = dofs.getElementBlock(block);
475  for(std::size_t e=0;e<elements.size();e++) {
476  std::vector<GO> gids;
477  dofs.getElementGIDs(elements[e],gids,block);
478 
479  // output gids belonging to this element
480  os << "[ " << gids[0];
481  for(std::size_t g=1;g<gids.size();g++)
482  os << ", " << gids[g];
483  os << " ]" << std::endl;
484  }
485  }
486 #endif
487 
488 }
489 
490 #endif
virtual int getNumFields() const =0
void getElementBlockIds(std::vector< std::string > &elementBlockIds) const
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Teuchos::RCP< const FieldPattern > getFieldPattern(const std::string &name) const
Find a field pattern stored for a particular block and field number. This will retrive the pattern ad...
virtual void getElementBlockIds(std::vector< std::string > &elementBlockIds) const =0
virtual const std::vector< int > & getBlockFieldNumbers(const std::string &blockId) const =0
PHX::MDField< const ScalarT, Cell, IP > field
const std::vector< Teuchos::RCP< UniqueGlobalIndexer< LocalOrdinalT, GlobalOrdinalT > > > & getFieldDOFManagers() const
virtual const std::string & getFieldString(int num) const =0
Reverse lookup of the field string from a field number.
Copy
bool fieldInBlock(const std::string &field, const std::string &block) const
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
#define TEUCHOS_ASSERT(assertion_test)
Teuchos::RCP< Thyra::LinearOpWithSolveFactoryBase< double > > buildLOWSFactory(bool blockedAssembly, const Teuchos::RCP< const panzer::UniqueGlobalIndexerBase > &globalIndexer, const Teuchos::RCP< panzer::ConnManagerBase< int > > &conn_manager, int spatialDim, const Teuchos::RCP< const Teuchos::MpiComm< int > > &mpi_comm, const Teuchos::RCP< Teuchos::ParameterList > &strat_params, bool writeCoordinates, bool writeTopo)
void writeTopology(const panzer::BlockedDOFManager< int, GO > &blkDofs)