1 // @HEADER 2 // ************************************************************************ 3 // 4 // Rapid Optimization Library (ROL) Package 5 // Copyright (2014) Sandia Corporation 6 // 7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive 8 // license for use of this work by or on behalf of the U.S. Government. 9 // 10 // Redistribution and use in source and binary forms, with or without 11 // modification, are permitted provided that the following conditions are 12 // met: 13 // 14 // 1. Redistributions of source code must retain the above copyright 15 // notice, this list of conditions and the following disclaimer. 16 // 17 // 2. Redistributions in binary form must reproduce the above copyright 18 // notice, this list of conditions and the following disclaimer in the 19 // documentation and/or other materials provided with the distribution. 20 // 21 // 3. Neither the name of the Corporation nor the names of the 22 // contributors may be used to endorse or promote products derived from 23 // this software without specific prior written permission. 24 // 25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY 26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR 28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE 29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 36 // 37 // Questions? Contact lead developers: 38 // Drew Kouri (dpkouri@sandia.gov) and 39 // Denis Ridzal (dridzal@sandia.gov) 40 // 41 // ************************************************************************ 42 // @HEADER 43 44 /*! \file filter.hpp 45 \brief Density filtering based on solving a PDE. 46 */ 47 48 #ifndef ROL_DENSITY_FILTER_H 49 #define ROL_DENSITY_FILTER_H 50 51 #include "Teuchos_GlobalMPISession.hpp" 52 #include "Teuchos_TimeMonitor.hpp" 53 54 #include "Tpetra_MultiVector.hpp" 55 #include "Tpetra_Vector.hpp" 56 #include "Tpetra_CrsGraph.hpp" 57 #include "Tpetra_CrsMatrix.hpp" 58 #include "Tpetra_Version.hpp" 59 #include "Tpetra_RowMatrixTransposer.hpp" 60 #include "MatrixMarket_Tpetra.hpp" 61 62 #include "Intrepid_HGRAD_QUAD_C2_FEM.hpp" 63 #include "Intrepid_HGRAD_QUAD_C1_FEM.hpp" 64 #include "Intrepid_DefaultCubatureFactory.hpp" 65 #include "Intrepid_FunctionSpaceTools.hpp" 66 #include "Intrepid_CellTools.hpp" 67 68 #include "Amesos2.hpp" 69 70 #include "../TOOLS/dofmanager.hpp" 71 72 template<class Real> 73 class DensityFilter { 74 75 private: 76 using GO = typename Tpetra::Map<>::global_ordinal_type; 77 78 ROL::Ptr<MeshManager<Real> > meshMgr_; 79 ROL::Ptr<DofManager<Real> > dofMgr_; 80 std::vector<ROL::Ptr<Intrepid::Basis<Real, Intrepid::FieldContainer<Real> > > > basisPtrs_; 81 82 ROL::Ptr<const Teuchos::Comm<int> > commPtr_; 83 int myRank_; 84 int numProcs_; 85 86 int basisOrder_; 87 88 ROL::Ptr<const Tpetra::Map<> > myOverlapMap_; 89 ROL::Ptr<const Tpetra::Map<> > myUniqueMap_; 90 ROL::Ptr<const Tpetra::Map<> > myBColumnMap_; 91 ROL::Ptr<Tpetra::CrsGraph<> > matAGraph_; 92 ROL::Ptr<Tpetra::CrsGraph<> > matBGraph_; 93 ROL::Ptr<Tpetra::CrsMatrix<> > matA_; 94 ROL::Ptr<Tpetra::CrsMatrix<> > matB_; 95 ROL::Ptr<Tpetra::CrsMatrix<> > matB_trans_; 96 ROL::Ptr<Tpetra::MultiVector<> > vecCellVolumes_; 97 98 Teuchos::Array<GO> myCellIds_; 99 Teuchos::Array<Real> myCellVolumes_; 100 101 ROL::Ptr<Amesos2::Solver< Tpetra::CrsMatrix<>, Tpetra::MultiVector<> > > solverA_; 102 103 shards::CellTopology cellType_; 104 int spaceDim_; 105 int numNodesPerCell_; 106 int numCubPoints_; 107 108 GO totalNumCells_; 109 GO totalNumDofs_; 110 GO numCells_; 111 112 ROL::Ptr<Intrepid::FieldContainer<Real> > cubPoints_; 113 ROL::Ptr<Intrepid::FieldContainer<Real> > cubWeights_; 114 ROL::Ptr<Intrepid::FieldContainer<Real> > cellNodes_; 115 ROL::Ptr<Intrepid::FieldContainer<Real> > cellJac_; 116 ROL::Ptr<Intrepid::FieldContainer<Real> > cellJacInv_; 117 ROL::Ptr<Intrepid::FieldContainer<Real> > cellJacDet_; 118 ROL::Ptr<Intrepid::FieldContainer<Real> > cellWeightedMeasure_; 119 ROL::Ptr<Intrepid::FieldContainer<Real> > valReference_; 120 ROL::Ptr<Intrepid::FieldContainer<Real> > gradReference_; 121 ROL::Ptr<Intrepid::FieldContainer<Real> > valPhysical_; 122 ROL::Ptr<Intrepid::FieldContainer<Real> > gradPhysical_; 123 ROL::Ptr<Intrepid::FieldContainer<Real> > kappaGradPhysical_; 124 ROL::Ptr<Intrepid::FieldContainer<Real> > valPhysicalWeighted_; 125 ROL::Ptr<Intrepid::FieldContainer<Real> > gradPhysicalWeighted_; 126 ROL::Ptr<Intrepid::FieldContainer<Real> > gradgradMats_; 127 ROL::Ptr<Intrepid::FieldContainer<Real> > valvalMats_; 128 ROL::Ptr<Intrepid::FieldContainer<Real> > valMats_; 129 ROL::Ptr<Intrepid::FieldContainer<Real> > onesVec_; 130 ROL::Ptr<Intrepid::FieldContainer<Real> > cubPointsPhysical_; 131 ROL::Ptr<Intrepid::FieldContainer<Real> > kappa_; 132 133 Real lengthScale_; 134 bool enableFilter_; 135 136 public: 137 DensityFilter(const ROL::Ptr<const Teuchos::Comm<int>> & comm,const Teuchos::RCP<Teuchos::ParameterList> & parlist,const ROL::Ptr<std::ostream> & outStream)138 DensityFilter(const ROL::Ptr<const Teuchos::Comm<int> > &comm, 139 const Teuchos::RCP<Teuchos::ParameterList> &parlist, 140 const ROL::Ptr<std::ostream> &outStream) { 141 142 lengthScale_ = parlist->sublist("Density Filter").get<Real>("Length Scale"); 143 lengthScale_ = std::pow(lengthScale_/static_cast<Real>(2*std::sqrt(3)), 2); 144 enableFilter_ = parlist->sublist("Density Filter").get<bool>("Enable"); 145 146 /************************************/ 147 /*** Retrieve communication data. ***/ 148 /************************************/ 149 commPtr_ = comm; 150 myRank_ = commPtr_->getRank(); 151 numProcs_ = commPtr_->getSize(); 152 *outStream << "Total number of processors: " << numProcs_ << std::endl; 153 /************************************/ 154 /************************************/ 155 156 /*************************************/ 157 /*** Retrieve parameter list data. ***/ 158 /*************************************/ 159 basisOrder_ = parlist->sublist("PDE FEM").get("Order of FE Discretization", 1); 160 int cellSplit = parlist->sublist("Geometry").get("Partition type", 1); 161 /*************************************/ 162 /*************************************/ 163 164 /****************************************************************************/ 165 /*** Initialize mesh / finite element fields / degree-of-freedom manager. ***/ 166 /****************************************************************************/ 167 168 // Mesh manager. 169 meshMgr_ = ROL::makePtr<MeshManager_Rectangle<Real>>(*parlist); 170 // Finite element fields. 171 ROL::Ptr<Intrepid::Basis<Real, Intrepid::FieldContainer<Real> > > basisPtr; 172 if (basisOrder_ == 1) { 173 basisPtr = ROL::makePtr<Intrepid::Basis_HGRAD_QUAD_C1_FEM<Real, Intrepid::FieldContainer<Real> >>(); 174 } 175 else if (basisOrder_ == 2) { 176 basisPtr = ROL::makePtr<Intrepid::Basis_HGRAD_QUAD_C2_FEM<Real, Intrepid::FieldContainer<Real> >>(); 177 } 178 basisPtrs_.resize(1, ROL::nullPtr); 179 basisPtrs_[0] = basisPtr; 180 // DOF coordinate interface. 181 ROL::Ptr<Intrepid::DofCoordsInterface<Intrepid::FieldContainer<Real> > > coord_iface = 182 ROL::dynamicPtrCast<Intrepid::DofCoordsInterface<Intrepid::FieldContainer<Real> > >(basisPtrs_[0]); 183 // Degree-of-freedom manager. 184 dofMgr_ = ROL::makePtr<DofManager<Real>>(meshMgr_, basisPtrs_); 185 // Retrieve total number of cells in the mesh. 186 totalNumCells_ = meshMgr_->getNumCells(); 187 // Retrieve total number of degrees of freedom in the mesh. 188 totalNumDofs_ = dofMgr_->getNumDofs(); 189 190 /****************************************************************************/ 191 /****************************************************************************/ 192 193 194 /****************************************************/ 195 /*** Build parallel communication infrastructure. ***/ 196 /****************************************************/ 197 198 // Partition the cells in the mesh. We use a basic quasi-equinumerous partitioning, 199 // where the remainder, if any, is assigned to the last processor. 200 Teuchos::Array<GO> myGlobIds_; 201 Teuchos::Array<GO> cellOffsets_(numProcs_, 0); 202 GO cellsPerProc = totalNumCells_ / numProcs_; 203 numCells_ = cellsPerProc; 204 switch(cellSplit) { 205 case 0: 206 if (myRank_ == 0) { // remainder in the first 207 numCells_ += totalNumCells_ % numProcs_; 208 } 209 for (int i=1; i<numProcs_; ++i) { 210 cellOffsets_[i] = cellOffsets_[i-1] + cellsPerProc + (static_cast<int>(i==1))*(totalNumCells_ % numProcs_); 211 } 212 break; 213 case 1: 214 if (myRank_ == numProcs_-1) { // remainder in the last 215 numCells_ += totalNumCells_ % numProcs_; 216 } 217 for (int i=1; i<numProcs_; ++i) { 218 cellOffsets_[i] = cellOffsets_[i-1] + cellsPerProc; 219 } 220 break; 221 case 2: 222 if (myRank_ < (totalNumCells_%numProcs_)) { // spread remainder, starting from the first 223 numCells_++; 224 } 225 for (int i=1; i<numProcs_; ++i) { 226 cellOffsets_[i] = cellOffsets_[i-1] + cellsPerProc + (static_cast<int>(i-1<(totalNumCells_%numProcs_))); 227 } 228 break; 229 } 230 Intrepid::FieldContainer<GO> &cellDofs = *(dofMgr_->getCellDofs()); 231 int numLocalDofs = cellDofs.dimension(1); 232 *outStream << "Cell offsets across processors: " << cellOffsets_ << std::endl; 233 for (GO i=0; i<numCells_; ++i) { 234 myCellIds_.push_back(cellOffsets_[myRank_]+i); 235 for (int j=0; j<numLocalDofs; ++j) { 236 myGlobIds_.push_back( cellDofs(cellOffsets_[myRank_]+i,j) ); 237 } 238 } 239 std::sort(myGlobIds_.begin(), myGlobIds_.end()); 240 myGlobIds_.erase( std::unique(myGlobIds_.begin(), myGlobIds_.end()), myGlobIds_.end() ); 241 242 // Build maps. 243 myOverlapMap_ = ROL::makePtr<Tpetra::Map<>>(Teuchos::OrdinalTraits<Tpetra::global_size_t>::invalid(), 244 myGlobIds_, 0, comm); 245 //std::cout << std::endl << myOverlapMap_->getNodeElementList(); 246 /** One can also use the non-member function: 247 myOverlapMap_ = Tpetra::createNonContigMap<int,int>(myGlobIds_, comm); 248 to build the overlap map. 249 **/ 250 myUniqueMap_ = Tpetra::createOneToOne(myOverlapMap_); 251 //std::cout << std::endl << myUniqueMap_->getNodeElementList() << std::endl; 252 253 myBColumnMap_ = ROL::makePtr<Tpetra::Map<>>(Teuchos::OrdinalTraits<Tpetra::global_size_t>::invalid(), 254 myCellIds_, 0, comm); 255 256 /****************************************************/ 257 /****************************************************/ 258 259 260 /****************************************************/ 261 /*** Set up local discretization data and arrays. ***/ 262 /****************************************************/ 263 264 // Retrieve some basic cell information. 265 cellType_ = (basisPtrs_[0])->getBaseCellTopology(); // get the cell type from any basis 266 spaceDim_ = cellType_.getDimension(); // retrieve spatial dimension 267 numNodesPerCell_ = cellType_.getNodeCount(); // retrieve number of nodes per cell 268 269 // Cubature data. 270 Intrepid::DefaultCubatureFactory<Real> cubFactory; // create cubature factory 271 int cubDegree = 4; // set cubature degree, e.g., 2 272 ROL::Ptr<Intrepid::Cubature<Real> > cellCub = cubFactory.create(cellType_, cubDegree); // create default cubature 273 numCubPoints_ = cellCub->getNumPoints(); // retrieve number of cubature points 274 275 int lfs = dofMgr_->getLocalFieldSize(0); 276 277 // Discretization data. 278 cubPoints_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCubPoints_, spaceDim_); 279 cubWeights_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCubPoints_); 280 cubPointsPhysical_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCells_, numCubPoints_, spaceDim_); 281 cellNodes_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCells_, numNodesPerCell_, spaceDim_); 282 cellJac_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCells_, numCubPoints_, spaceDim_, spaceDim_); 283 cellJacInv_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCells_, numCubPoints_, spaceDim_, spaceDim_); 284 cellJacDet_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCells_, numCubPoints_); 285 cellWeightedMeasure_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCells_, numCubPoints_); 286 valReference_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(lfs, numCubPoints_); 287 gradReference_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(lfs, numCubPoints_, spaceDim_); 288 valPhysical_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCells_, lfs, numCubPoints_); 289 gradPhysical_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCells_, lfs, numCubPoints_, spaceDim_); 290 kappaGradPhysical_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCells_, lfs, numCubPoints_, spaceDim_); 291 valPhysicalWeighted_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCells_, lfs, numCubPoints_); 292 gradPhysicalWeighted_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCells_, lfs, numCubPoints_, spaceDim_); 293 gradgradMats_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCells_, lfs, lfs); 294 valvalMats_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCells_, lfs, lfs); 295 valMats_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCells_, lfs, 1); 296 onesVec_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCells_, 1, numCubPoints_); 297 kappa_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCells_, numCubPoints_); 298 299 // Geometric definition of the cells in the mesh, based on the cell-to-node map and the domain partition. 300 Intrepid::FieldContainer<Real> &nodes = *meshMgr_->getNodes(); 301 Intrepid::FieldContainer<GO> &ctn = *meshMgr_->getCellToNodeMap(); 302 for (GO i=0; i<numCells_; ++i) { 303 for (int j=0; j<numNodesPerCell_; ++j) { 304 for (int k=0; k<spaceDim_; ++k) { 305 (*cellNodes_)(i, j, k) = nodes(ctn(myCellIds_[i],j), k); 306 } 307 } 308 } 309 310 /****************************************************/ 311 /****************************************************/ 312 313 314 /****************************************************************/ 315 /*** Assemble cellwise contributions to vectors and matrices. ***/ 316 /****************************************************************/ 317 318 cellCub->getCubature(*cubPoints_, *cubWeights_); // retrieve cubature points and weights 319 (*basisPtrs_[0]).getValues(*gradReference_, *cubPoints_, Intrepid::OPERATOR_GRAD); // evaluate grad operator at cubature points 320 (*basisPtrs_[0]).getValues(*valReference_, *cubPoints_, Intrepid::OPERATOR_VALUE); // evaluate value operator at cubature points 321 322 Intrepid::CellTools<Real>::setJacobian(*cellJac_, *cubPoints_, *cellNodes_, cellType_); // compute cell Jacobians 323 Intrepid::CellTools<Real>::setJacobianInv(*cellJacInv_, *cellJac_); // compute inverses of cell Jacobians 324 Intrepid::CellTools<Real>::setJacobianDet(*cellJacDet_, *cellJac_); // compute determinants of cell Jacobians 325 326 Intrepid::FunctionSpaceTools::computeCellMeasure<Real>(*cellWeightedMeasure_, // compute weighted cell measure 327 *cellJacDet_, 328 *cubWeights_); 329 330 Intrepid::CellTools<Real>::mapToPhysicalFrame(*cubPointsPhysical_, // map reference cubature points to physical space 331 *cubPoints_, 332 *cellNodes_, 333 cellType_); 334 335 Intrepid::FunctionSpaceTools::HGRADtransformGRAD<Real>(*gradPhysical_, // transform reference gradients into physical space 336 *cellJacInv_, 337 *gradReference_); 338 Intrepid::FunctionSpaceTools::multiplyMeasure<Real>(*gradPhysicalWeighted_, // multiply with weighted measure 339 *cellWeightedMeasure_, 340 *gradPhysical_); 341 for (int i=0; i<numCells_; ++i) { // evaluate conductivity kappa at cubature points 342 for (int j=0; j<numCubPoints_; ++j) { 343 (*kappa_)(i, j) = funcKappa((*cubPointsPhysical_)(i, j, 0), 344 (*cubPointsPhysical_)(i, j, 1)); 345 } 346 } 347 Intrepid::FunctionSpaceTools::tensorMultiplyDataField<Real>(*kappaGradPhysical_, // multiply with conductivity kappa 348 *kappa_, 349 *gradPhysical_); 350 Intrepid::FunctionSpaceTools::integrate<Real>(*gradgradMats_, // compute local grad.(kappa)grad (stiffness) matrices 351 *kappaGradPhysical_, 352 *gradPhysicalWeighted_, 353 Intrepid::COMP_CPP); 354 355 Intrepid::FunctionSpaceTools::HGRADtransformVALUE<Real>(*valPhysical_, // transform reference values into physical space 356 *valReference_); 357 Intrepid::FunctionSpaceTools::multiplyMeasure<Real>(*valPhysicalWeighted_, // multiply with weighted measure 358 *cellWeightedMeasure_, 359 *valPhysical_); 360 Intrepid::FunctionSpaceTools::integrate<Real>(*valvalMats_, // compute local val.val (mass) matrices 361 *valPhysical_, 362 *valPhysicalWeighted_, 363 Intrepid::COMP_CPP); 364 365 for (int i=0; i<numCells_; ++i) { // fill vector of ones 366 for (int j=0; j<numCubPoints_; ++j) { 367 (*onesVec_)(i, 0, j) = static_cast<Real>(1); 368 } 369 } 370 Intrepid::FunctionSpaceTools::integrate<Real>(*valMats_, // compute local val.1 matrices 371 *valPhysicalWeighted_, 372 *onesVec_, 373 Intrepid::COMP_CPP); 374 375 /****************************************************************/ 376 /****************************************************************/ 377 378 /****************************************/ 379 /*** Assemble global data structures. ***/ 380 /****************************************/ 381 382 // Assemble graphs. 383 matAGraph_ = ROL::makePtr<Tpetra::CrsGraph<>>(myUniqueMap_, 0); 384 Teuchos::ArrayRCP<const GO> cellDofsArrayRCP = cellDofs.getData(); 385 for (int i=0; i<numCells_; ++i) { 386 for (int j=0; j<numLocalDofs; ++j) { 387 matAGraph_->insertGlobalIndices(cellDofs(myCellIds_[i],j), cellDofsArrayRCP(myCellIds_[i]*numLocalDofs, numLocalDofs)); 388 } 389 } 390 matAGraph_->fillComplete(); 391 matBGraph_ = ROL::makePtr<Tpetra::CrsGraph<>>(myUniqueMap_, myBColumnMap_, 0); 392 Teuchos::ArrayRCP<const GO> cellIdsArrayRCP = Teuchos::arcpFromArray(myCellIds_); 393 for (int i=0; i<numCells_; ++i) { 394 for (int j=0; j<numLocalDofs; ++j) { 395 matBGraph_->insertGlobalIndices(cellDofs(myCellIds_[i],j), cellIdsArrayRCP(i, 1)); 396 } 397 } 398 matBGraph_->fillComplete(myBColumnMap_, myUniqueMap_); 399 400 // Assemble matrices. 401 // Filter matrix = stiffness matrix plus mass matrix. 402 matA_ = ROL::makePtr<Tpetra::CrsMatrix<>>(matAGraph_); 403 int numLocalMatEntries = numLocalDofs*numLocalDofs; 404 Teuchos::ArrayRCP<const Real> gradgradArrayRCP = gradgradMats_->getData(); 405 Teuchos::ArrayRCP<const Real> valvalArrayRCP = valvalMats_->getData(); 406 for (int i=0; i<numCells_; ++i) { 407 for (int j=0; j<numLocalDofs; ++j) { 408 matA_->sumIntoGlobalValues(cellDofs(myCellIds_[i],j), 409 cellDofsArrayRCP(myCellIds_[i]*numLocalDofs, numLocalDofs), 410 gradgradArrayRCP(i*numLocalMatEntries+j*numLocalDofs, numLocalDofs)); 411 matA_->sumIntoGlobalValues(cellDofs(myCellIds_[i],j), 412 cellDofsArrayRCP(myCellIds_[i]*numLocalDofs, numLocalDofs), 413 valvalArrayRCP(i*numLocalMatEntries+j*numLocalDofs, numLocalDofs)); 414 } 415 } 416 matA_->fillComplete(); 417 // B matrix. 418 matB_ = ROL::makePtr<Tpetra::CrsMatrix<>>(matBGraph_); 419 Teuchos::ArrayRCP<const Real> valArrayRCP = valMats_->getData(); 420 for (int i=0; i<numCells_; ++i) { 421 for (int j=0; j<numLocalDofs; ++j) { 422 matB_->sumIntoGlobalValues(cellDofs(myCellIds_[i],j), 423 cellIdsArrayRCP(i, 1), 424 valArrayRCP(i*numLocalDofs+j, 1)); 425 } 426 } 427 matB_->fillComplete(); 428 429 // Compute cell colume vector. 430 computeCellVolumes(); 431 432 // Create matrix transposes. 433 Tpetra::RowMatrixTransposer<> transposerB(matB_); 434 matB_trans_ = transposerB.createTranspose(); 435 436 /*********************************/ 437 /*** Construct solver objects. ***/ 438 /*********************************/ 439 440 // Construct solver using Amesos2 factory. 441 try{ 442 solverA_ = Amesos2::create< Tpetra::CrsMatrix<>,Tpetra::MultiVector<> >("KLU2", matA_); 443 } catch (std::invalid_argument& e) { 444 std::cout << e.what() << std::endl; 445 } 446 solverA_->numericFactorization(); 447 448 /****************************************/ 449 /****************************************/ 450 451 //outputTpetraData(); 452 453 } 454 apply(ROL::Ptr<Tpetra::MultiVector<>> & Fx,const ROL::Ptr<const Tpetra::MultiVector<>> & x)455 void apply(ROL::Ptr<Tpetra::MultiVector<> > & Fx, const ROL::Ptr<const Tpetra::MultiVector<> > & x) { 456 if (enableFilter_) { 457 ROL::Ptr<Tpetra::MultiVector<> > Bx = ROL::makePtr<Tpetra::MultiVector<>>(matB_->getRangeMap(), 1); 458 ROL::Ptr<Tpetra::MultiVector<> > AiBx = ROL::makePtr<Tpetra::MultiVector<>>(matA_->getDomainMap(), 1); 459 ROL::Ptr<Tpetra::MultiVector<> > Fx_unscaled = ROL::makePtr<Tpetra::MultiVector<>>(matB_trans_->getRangeMap(), 1); 460 matB_->apply(*x, *Bx); 461 solverA_->setX(AiBx); 462 solverA_->setB(Bx); 463 solverA_->solve(); 464 //outputTpetraVector(AiBx, "density_nodal.txt"); 465 matB_trans_->apply(*AiBx, *Fx_unscaled); 466 ROL::Ptr<Tpetra::MultiVector<> > vecInvCellVolumes = ROL::makePtr<Tpetra::MultiVector<>>(matB_->getDomainMap(), 1); 467 vecInvCellVolumes->reciprocal(*vecCellVolumes_); 468 Fx->elementWiseMultiply(1.0, *(vecInvCellVolumes->getVector(0)), *Fx_unscaled, 0.0); 469 } 470 else { 471 Fx->update(1.0, *x, 0.0); 472 } 473 } 474 getMatA() const475 ROL::Ptr<Tpetra::CrsMatrix<> > getMatA() const { 476 return matA_; 477 } 478 getMatB(const bool & transpose=false) const479 ROL::Ptr<Tpetra::CrsMatrix<> > getMatB(const bool &transpose = false) const { 480 if (transpose) { 481 return matB_trans_; 482 } 483 else { 484 return matB_; 485 } 486 } 487 getSolver() const488 ROL::Ptr<Amesos2::Solver< Tpetra::CrsMatrix<>, Tpetra::MultiVector<> > > getSolver() const { 489 return solverA_; 490 } 491 funcKappa(const Real & x1,const Real & x2) const492 Real funcKappa(const Real &x1, const Real &x2) const { 493 return lengthScale_; 494 } 495 computeCellVolumes()496 void computeCellVolumes() { 497 for (int i=0; i<numCells_; ++i){ 498 Real temp = 0.0; 499 for(int j=0; j<numCubPoints_; ++j){ 500 temp += (*cellWeightedMeasure_)(i, j); 501 } 502 myCellVolumes_.push_back(temp); 503 } 504 505 vecCellVolumes_ = ROL::makePtr<Tpetra::MultiVector<>>(matB_->getDomainMap(), 1, true); 506 for (int i=0; i<numCells_; ++i){ 507 vecCellVolumes_->replaceGlobalValue(myCellIds_[i], 0, myCellVolumes_[i]); 508 } 509 } 510 outputTpetraData() const511 void outputTpetraData() const { 512 Tpetra::MatrixMarket::Writer< Tpetra::CrsMatrix<> > matWriter; 513 matWriter.writeSparseFile("filter_mat", matA_, true); 514 matWriter.writeSparseFile("projection_mat", matB_, true); 515 matWriter.writeSparseFile("projection_trans_mat", matB_trans_, true); 516 } 517 outputTpetraVector(const ROL::Ptr<const Tpetra::MultiVector<>> & vec,const std::string & filename) const518 void outputTpetraVector(const ROL::Ptr<const Tpetra::MultiVector<> > &vec, 519 const std::string &filename) const { 520 Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<> > vecWriter; 521 vecWriter.writeDenseFile(filename, vec); 522 } 523 524 }; // class StefanBoltzmannData 525 526 #endif 527