Albany: a Trilinos-based PDE code: Main_MeshSample.cpp Source File

Go to the documentation of this file.
00001 //*****************************************************************//
00002 //    Albany 2.0:  Copyright 2012 Sandia Corporation               //
00003 //    This Software is released under the BSD license detailed     //
00004 //    in the file "license.txt" in the top-level Albany directory  //
00005 //*****************************************************************//
00006 
00007 #include "Albany_DataTypes.hpp"
00008 #include "Albany_Utils.hpp"
00009 
00010 #include "Albany_DiscretizationFactory.hpp"
00011 #include "Albany_STKDiscretization.hpp"
00012 #include "Albany_AbstractSTKMeshStruct.hpp"
00013 
00014 #include "Albany_MORDiscretizationUtils.hpp"
00015 #include "Albany_StkNodalBasisSource.hpp"
00016 
00017 #include "MOR_WindowedAtomicBasisSource.hpp"
00018 #include "MOR_GreedyAtomicBasisSample.hpp"
00019 #include "MOR_StkNodalMeshReduction.hpp"
00020 #include "MOR_CollocationMetricCriterionFactory.hpp"
00021 
00022 #include "Epetra_Comm.h"
00023 #include "Epetra_Vector.h"
00024 #include "Epetra_Import.h"
00025 
00026 #include "Teuchos_Ptr.hpp"
00027 #include "Teuchos_RCP.hpp"
00028 #include "Teuchos_ArrayRCP.hpp"
00029 #include "Teuchos_Array.hpp"
00030 #include "Teuchos_Tuple.hpp"
00031 #include "Teuchos_GlobalMPISession.hpp"
00032 #include "Teuchos_Comm.hpp"
00033 #include "Teuchos_FancyOStream.hpp"
00034 #include "Teuchos_VerboseObject.hpp"
00035 #include "Teuchos_ParameterList.hpp"
00036 #include "Teuchos_XMLParameterListHelpers.hpp"
00037 #include "Teuchos_TestForException.hpp"
00038 
00039 #include <string>
00040 
00041 using Teuchos::RCP;
00042 
00043 class SampleDiscretization : public Albany::DiscretizationTransformation {
00044 public:
00045   SampleDiscretization(
00046       const Teuchos::ArrayView<const stk::mesh::EntityId> &nodeIds,
00047       const Teuchos::ArrayView<const stk::mesh::EntityId> &sensorNodeIds,
00048       bool performReduction);
00049 
00050   void operator()(Albany::DiscretizationFactory &discFactory);
00051 
00052 private:
00053   Teuchos::ArrayView<const stk::mesh::EntityId> nodeIds_;
00054   Teuchos::ArrayView<const stk::mesh::EntityId> sensorNodeIds_;
00055   bool performReduction_;
00056 };
00057 
00058 SampleDiscretization::SampleDiscretization(
00059     const Teuchos::ArrayView<const stk::mesh::EntityId> &nodeIds,
00060     const Teuchos::ArrayView<const stk::mesh::EntityId> &sensorNodeIds,
00061     bool performReduction) :
00062   nodeIds_(nodeIds),
00063   sensorNodeIds_(sensorNodeIds),
00064   performReduction_(performReduction)
00065 {}
00066 
00067 void
00068 SampleDiscretization::operator()(Albany::DiscretizationFactory &discFactory)
00069 {
00070   const RCP<Albany::AbstractMeshStruct> meshStruct = discFactory.getMeshStruct();
00071   const RCP<Albany::AbstractSTKMeshStruct> stkMeshStruct =
00072     Teuchos::rcp_dynamic_cast<Albany::AbstractSTKMeshStruct>(meshStruct, /*throw_on_fail =*/ true);
00073 
00074   stk::mesh::BulkData &bulkData = *stkMeshStruct->bulkData;
00075 
00076   stk::mesh::Part &samplePart = *stkMeshStruct->nsPartVec["sample_nodes"];
00077   MOR::addNodesToPart(nodeIds_, samplePart, bulkData);
00078 
00079   if (sensorNodeIds_.size() > 0) {
00080     stk::mesh::Part &sensorPart = *stkMeshStruct->nsPartVec["sensors"];
00081     MOR::addNodesToPart(sensorNodeIds_, sensorPart, bulkData);
00082   }
00083 
00084   if (performReduction_) {
00085     MOR::performNodalMeshReduction(samplePart, bulkData);
00086   }
00087 }
00088 
00089 RCP<Albany::AbstractDiscretization> sampledDiscretizationNew(
00090     const RCP<Teuchos::ParameterList> &topLevelParams,
00091     const RCP<const Epetra_Comm> &epetraComm,
00092     const Teuchos::ArrayView<const stk::mesh::EntityId> &nodeIds,
00093     const Teuchos::ArrayView<const stk::mesh::EntityId> &sensorNodeIds,
00094     bool performReduction)
00095 {
00096   SampleDiscretization transformation(nodeIds, sensorNodeIds, performReduction);
00097   return Albany::modifiedDiscretizationNew(topLevelParams, epetraComm, transformation);
00098 }
00099 
00100 void transferSolutionHistoryImpl(
00101     Albany::STKDiscretization &source,
00102     Albany::AbstractDiscretization &target,
00103     int depth)
00104 {
00105   Epetra_Vector targetVec(*target.getOverlapMap(), false);
00106   Epetra_Import importer(targetVec.Map(), *source.getMap());
00107 
00108   const RCP<Albany::AbstractSTKMeshStruct> sourceMeshStruct = source.getSTKMeshStruct();
00109 
00110   for (int rank = 0; rank != depth; ++rank) {
00111     const double stamp = sourceMeshStruct->getSolutionFieldHistoryStamp(rank);
00112     sourceMeshStruct->loadSolutionFieldHistory(rank);
00113     const RCP<const Epetra_Vector> sourceVec = source.getSolutionField();
00114     targetVec.Import(*sourceVec, importer, Insert);
00115     target.writeSolution(targetVec, stamp, /*overlapped =*/ true);
00116   }
00117 }
00118 
00119 void transferSolutionHistory(
00120     Albany::STKDiscretization &source,
00121     Albany::AbstractDiscretization &target)
00122 {
00123   const RCP<Albany::AbstractSTKMeshStruct> sourceMeshStruct = source.getSTKMeshStruct();
00124   const int steps = sourceMeshStruct->getSolutionFieldHistoryDepth();
00125   transferSolutionHistoryImpl(source, target, steps);
00126 }
00127 
00128 void transferSolutionHistory(
00129     Albany::STKDiscretization &source,
00130     Albany::AbstractDiscretization &target,
00131     int depthMax)
00132 {
00133   const RCP<Albany::AbstractSTKMeshStruct> sourceMeshStruct = source.getSTKMeshStruct();
00134   const int steps = std::min(sourceMeshStruct->getSolutionFieldHistoryDepth(), depthMax);
00135   transferSolutionHistoryImpl(source, target, steps);
00136 }
00137 
00138 RCP<Teuchos::ParameterEntry> getEntryCopy(
00139     const Teuchos::ParameterList &l,
00140     const std::string &name)
00141 {
00142   const Teuchos::Ptr<const Teuchos::ParameterEntry> source(l.getEntryPtr(name));
00143   if (Teuchos::nonnull(source)) {
00144     return Teuchos::rcp(new Teuchos::ParameterEntry(*source));
00145   }
00146   return Teuchos::null;
00147 }
00148 
00149 int main(int argc, char *argv[])
00150 {
00151   // Communicators
00152   Teuchos::GlobalMPISession mpiSession(&argc, &argv);
00153   const Albany_MPI_Comm nativeComm = Albany_MPI_COMM_WORLD;
00154   const RCP<const Teuchos::Comm<int> > teuchosComm = Albany::createTeuchosCommFromMpiComm(nativeComm);
00155   const RCP<const Epetra_Comm> epetraComm = Albany::createEpetraCommFromMpiComm(nativeComm);
00156 
00157   // Standard output
00158   const RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
00159 
00160   // Parse command-line argument for input file
00161   const std::string firstArg = (argc > 1) ? argv[1] : "";
00162   if (firstArg.empty() || firstArg == "--help") {
00163     *out << "AlbanyRBGen input-file-path\n";
00164     return 0;
00165   }
00166   const std::string inputFileName = argv[1];
00167 
00168   // Parse XML input file
00169   const RCP<Teuchos::ParameterList> topLevelParams = Teuchos::createParameterList("Albany Parameters");
00170   Teuchos::updateParametersFromXmlFileAndBroadcast(inputFileName, topLevelParams.ptr(), *teuchosComm);
00171 
00172   const bool sublistMustExist = true;
00173 
00174   // Setup discretization factory
00175   const RCP<Teuchos::ParameterList> discParams = Teuchos::sublist(topLevelParams, "Discretization", sublistMustExist);
00176   TEUCHOS_TEST_FOR_EXCEPT(discParams->get<std::string>("Method") != "Ioss");
00177   const std::string outputParamLabel = "Exodus Output File Name";
00178   const std::string sampledOutputParamLabel = "Reference Exodus Output File Name";
00179   const RCP<const Teuchos::ParameterEntry> reducedOutputParamEntry = getEntryCopy(*discParams, outputParamLabel);
00180   const RCP<const Teuchos::ParameterEntry> sampledOutputParamEntry = getEntryCopy(*discParams, sampledOutputParamLabel);
00181   discParams->remove(outputParamLabel, /*throwIfNotExists =*/ false);
00182   discParams->remove(sampledOutputParamLabel, /*throwIfNotExists =*/ false);
00183   const RCP<const Teuchos::ParameterList> discParamsCopy = Teuchos::rcp(new Teuchos::ParameterList(*discParams));
00184 
00185   const RCP<Teuchos::ParameterList> problemParams = Teuchos::sublist(topLevelParams, "Problem", sublistMustExist);
00186   const RCP<const Teuchos::ParameterList> problemParamsCopy = Teuchos::rcp(new Teuchos::ParameterList(*problemParams));
00187 
00188   // Create original (full) discretization
00189   const RCP<Albany::AbstractDiscretization> disc = Albany::discretizationNew(topLevelParams, epetraComm);
00190 
00191   // Determine mesh sample
00192   const RCP<Teuchos::ParameterList> samplingParams = Teuchos::sublist(topLevelParams, "Mesh Sampling", sublistMustExist);
00193   const int firstVectorRank = samplingParams->get("First Vector Rank", 0);
00194   const Teuchos::Ptr<const int> basisSizeMax = Teuchos::ptr(samplingParams->getPtr<int>("Basis Size Max"));
00195   const int sampleSize = samplingParams->get("Sample Size", 0);
00196 
00197   *out << "Sampling " << sampleSize << " nodes";
00198   if (Teuchos::nonnull(basisSizeMax)) {
00199     *out << " based on no more than " << *basisSizeMax << " basis vectors";
00200   }
00201   if (firstVectorRank != 0) {
00202     *out << " starting from vector rank " << firstVectorRank;
00203   }
00204   *out << "\n";
00205 
00206   const RCP<Albany::STKDiscretization> stkDisc =
00207     Teuchos::rcp_dynamic_cast<Albany::STKDiscretization>(disc, /*throw_on_fail =*/ true);
00208   const RCP<MOR::AtomicBasisSource> rawBasisSource = Teuchos::rcp(new Albany::StkNodalBasisSource(stkDisc));
00209   const RCP<MOR::AtomicBasisSource> basisSource = Teuchos::rcp(
00210       Teuchos::nonnull(basisSizeMax) ?
00211       new MOR::WindowedAtomicBasisSource(rawBasisSource, firstVectorRank, *basisSizeMax) :
00212       new MOR::WindowedAtomicBasisSource(rawBasisSource, firstVectorRank)
00213       );
00214 
00215   MOR::CollocationMetricCriterionFactory criterionFactory(samplingParams);
00216   const Teuchos::RCP<const MOR::CollocationMetricCriterion> criterion =
00217     criterionFactory.instanceNew(basisSource->entryCountMax());
00218   const Teuchos::RCP<MOR::GreedyAtomicBasisSample> sampler(new MOR::GreedyAtomicBasisSample(*basisSource, criterion));
00219   sampler->sampleSizeInc(sampleSize);
00220 
00221   Teuchos::Array<stk::mesh::EntityId> sampleNodeIds;
00222   const Teuchos::ArrayView<const int> sampleAtoms = sampler->sample();
00223   sampleNodeIds.reserve(sampleAtoms.size());
00224   for (Teuchos::ArrayView<const int>::const_iterator it = sampleAtoms.begin(), it_end = sampleAtoms.end(); it != it_end; ++it) {
00225     sampleNodeIds.push_back(*it + 1);
00226   }
00227 
00228   *out << "Sample = " << sampleNodeIds << "\n";
00229 
00230   // Choose first sample node as sensor
00231   const Teuchos::ArrayView<const stk::mesh::EntityId> sensorNodeIds = sampleNodeIds.view(0, 1);
00232 
00233   const Teuchos::Array<std::string> additionalNodeSets =
00234     Teuchos::tuple(std::string("sample_nodes"), std::string("sensors"));
00235 
00236   // Create sampled discretization
00237   if (Teuchos::nonnull(sampledOutputParamEntry)) {
00238     const RCP<Teuchos::ParameterList> discParamsLocalCopy = Teuchos::rcp(new Teuchos::ParameterList(*discParamsCopy));
00239     discParamsLocalCopy->setEntry("Exodus Output File Name", *sampledOutputParamEntry);
00240     discParamsLocalCopy->set("Additional Node Sets", additionalNodeSets);
00241     topLevelParams->set("Discretization", *discParamsLocalCopy);
00242     topLevelParams->set("Problem", *problemParamsCopy);
00243 
00244     const bool performReduction = false;
00245     const RCP<Albany::AbstractDiscretization> sampledDisc =
00246       sampledDiscretizationNew(topLevelParams, epetraComm, sampleNodeIds, sensorNodeIds, performReduction);
00247 
00248     if (Teuchos::nonnull(basisSizeMax)) {
00249       transferSolutionHistory(*stkDisc, *sampledDisc, *basisSizeMax + firstVectorRank);
00250     } else {
00251       transferSolutionHistory(*stkDisc, *sampledDisc);
00252     }
00253   }
00254 
00255   // Create reduced discretization
00256   if (Teuchos::nonnull(reducedOutputParamEntry)) {
00257     const RCP<Teuchos::ParameterList> discParamsLocalCopy = Teuchos::rcp(new Teuchos::ParameterList(*discParamsCopy));
00258     discParamsLocalCopy->setEntry("Exodus Output File Name", *reducedOutputParamEntry);
00259     discParamsLocalCopy->set("Additional Node Sets", additionalNodeSets);
00260     topLevelParams->set("Discretization", *discParamsLocalCopy);
00261     topLevelParams->set("Problem", *problemParamsCopy);
00262 
00263     const bool performReduction = true;
00264     const RCP<Albany::AbstractDiscretization> reducedDisc =
00265       sampledDiscretizationNew(topLevelParams, epetraComm, sampleNodeIds, sensorNodeIds, performReduction);
00266 
00267     if (Teuchos::nonnull(basisSizeMax)) {
00268       transferSolutionHistory(*stkDisc, *reducedDisc, *basisSizeMax + firstVectorRank);
00269     } else {
00270       transferSolutionHistory(*stkDisc, *reducedDisc);
00271     }
00272   }
00273 }