Albany: a Trilinos-based PDE code: QCAD_ResponseFieldIntegral

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 <fstream>
00008 #include "Teuchos_TestForException.hpp"
00009 #include "Teuchos_CommHelpers.hpp"
00010 #include "Albany_Utils.hpp"
00011 
00012 template<typename EvalT, typename Traits>
00013 QCAD::ResponseFieldIntegral<EvalT, Traits>::
00014 ResponseFieldIntegral(Teuchos::ParameterList& p,
00015           const Teuchos::RCP<Albany::Layouts>& dl) :
00016   coordVec("Coord Vec", dl->qp_vector),
00017   weights("Weights", dl->qp_scalar)
00018 {
00020   Teuchos::ParameterList* plist = 
00021     p.get<Teuchos::ParameterList*>("Parameter List");
00022   Teuchos::RCP<const Teuchos::ParameterList> reflist = 
00023     this->getValidResponseParameters();
00024   plist->validateParameters(*reflist,0);
00025 
00027   Teuchos::RCP<Teuchos::ParameterList> paramsFromProblem = 
00028     p.get< Teuchos::RCP<Teuchos::ParameterList> >("Parameters From Problem");
00029   if(paramsFromProblem != Teuchos::null) {
00030 
00031     // Material database 
00032     materialDB = paramsFromProblem->get< Teuchos::RCP<QCAD::MaterialDatabase> >("MaterialDB");
00033 
00034     // Length unit in meters
00035     length_unit_in_m = paramsFromProblem->get<double>("Length unit in m");
00036   }
00037   else {
00038     materialDB = Teuchos::null;
00039     length_unit_in_m = 1.0e-6; //default length unit = microns (backward compat)
00040   }
00041        
00043   Teuchos::RCP<PHX::DataLayout> scalar_dl = dl->qp_scalar;
00044   numQPs = scalar_dl->dimension(1);
00045   
00047   Teuchos::RCP<PHX::DataLayout> vector_dl = dl->qp_vector;
00048   std::vector<PHX::DataLayout::size_type> dims;
00049   vector_dl->dimensions(dims);
00050   numDims = dims[2];
00051 
00053   opRegion  = Teuchos::rcp( new QCAD::MeshRegion<EvalT, Traits>("Coord Vec","Weights",*plist,materialDB,dl) );
00054 
00056   std::string fieldName;
00057 
00058   fieldName = plist->get<std::string>("Field Name","");
00059   if(fieldName.length() > 0) {
00060     fieldNames.push_back(fieldName);
00061     conjugateFieldFlag.push_back(plist->get<bool>("Conjugate Field",false));
00062 
00063     fieldName = plist->get<std::string>("Field Name Im","");
00064     fieldNames_Imag.push_back(fieldName);
00065     if(fieldName.length() > 0) fieldIsComplex.push_back(true);
00066     else fieldIsComplex.push_back(false);
00067   }
00068 
00069   for(int i=1; i < QCAD::MAX_FIELDNAMES_IN_INTEGRAL; i++) {
00070     fieldName = plist->get<std::string>(Albany::strint("Field Name",i),"");
00071     if(fieldName.length() > 0) {
00072       fieldNames.push_back(fieldName);
00073       conjugateFieldFlag.push_back(plist->get<bool>(Albany::strint("Conjugate Field",i),false));
00074 
00075       fieldName = plist->get<std::string>(Albany::strint("Field Name Im",i),"");
00076       fieldNames_Imag.push_back(fieldName);
00077       if(fieldName.length() > 0) fieldIsComplex.push_back(true);
00078       else fieldIsComplex.push_back(false);
00079     }
00080     else break;
00081   }
00082   bReturnImagPart = plist->get<bool>("Return Imaginary Part",false);
00083   
00084   std::string integrandLinLengthUnit; // linear length unit of integrand (e.g. "cm" for integrand in cm^-3)
00085   integrandLinLengthUnit = plist->get<std::string>("Integrand Length Unit","cm");
00086   bPositiveOnly = plist->get<bool>("Positive Return Only",false);
00087 
00089   double integrand_length_unit_in_m;
00090   if( integrandLinLengthUnit == "m" )       integrand_length_unit_in_m = 1.0;
00091   else if( integrandLinLengthUnit == "cm" ) integrand_length_unit_in_m = 1e-2;
00092   else if( integrandLinLengthUnit == "um" ) integrand_length_unit_in_m = 1e-6;
00093   else if( integrandLinLengthUnit == "nm" ) integrand_length_unit_in_m = 1e-9;
00094   else if( integrandLinLengthUnit == "mesh" ) integrand_length_unit_in_m = length_unit_in_m;
00095   else integrand_length_unit_in_m = length_unit_in_m;  // assume same unit as mesh (e.g. if unit string is blank)
00096   
00097   double X0 = length_unit_in_m / integrand_length_unit_in_m; // length scaling to get to integrand's lenght unit
00098 
00099   if (numDims == 1)       scaling = X0; 
00100   else if (numDims == 2)  scaling = X0*X0; 
00101   else if (numDims == 3)  scaling = X0*X0*X0; 
00102   else 
00103     TEUCHOS_TEST_FOR_EXCEPTION (true, Teuchos::Exceptions::InvalidParameter, std::endl 
00104         << "Error! Invalid number of dimensions: " << numDims << std::endl);
00105 
00106 
00108   std::vector<std::string>::const_iterator it;
00109   //for(it = fieldNames.begin(); it != fieldNames.end(); ++it) {
00110   for(std::size_t i=0; i<fieldNames.size(); i++) {
00111     PHX::MDField<ScalarT,Cell,QuadPoint> f(fieldNames[i], scalar_dl);
00112     fields.push_back(f); this->addDependentField(f);
00113 
00114     PHX::MDField<ScalarT,Cell,QuadPoint> fi(fieldNames_Imag[i], scalar_dl);
00115     fields_Imag.push_back(fi);
00116 
00117     if(fieldIsComplex[i]) this->addDependentField(fi);
00118   }
00119 
00120   this->addDependentField(coordVec);
00121   this->addDependentField(weights);
00122   opRegion->addDependentFields(this);
00123 
00124   //TODO: make name unique? Is this needed for anything?
00125   this->setName(fieldName+" Response Field Integral"+PHX::TypeString<EvalT>::value);
00126   
00127   // Setup scatter evaluator
00128   p.set("Stand-alone Evaluator", false);
00129   std::string local_response_name = 
00130     fieldName + " Local Response Field Integral";
00131   std::string global_response_name = 
00132     fieldName + " Global Response Field Integral";
00133   PHX::Tag<ScalarT> local_response_tag(local_response_name, 
00134                dl->cell_scalar);
00135   PHX::Tag<ScalarT> global_response_tag(global_response_name, 
00136           dl->workset_scalar);
00137   p.set("Local Response Field Tag", local_response_tag);
00138   p.set("Global Response Field Tag", global_response_tag);
00139   PHAL::SeparableScatterScalarResponse<EvalT,Traits>::setup(p,dl);
00140 }
00141 
00142 // **********************************************************************
00143 template<typename EvalT, typename Traits>
00144 void QCAD::ResponseFieldIntegral<EvalT, Traits>::
00145 postRegistrationSetup(typename Traits::SetupData d,
00146                       PHX::FieldManager<Traits>& fm)
00147 {
00148   typename std::vector<PHX::MDField<ScalarT,Cell,QuadPoint> >::iterator it;
00149   //for(it = fields.begin(); it != fields.end(); ++it)
00150   for(std::size_t i=0; i<fields.size(); i++) {
00151     this->utils.setFieldData(fields[i],fm);
00152     if(fieldIsComplex[i]) this->utils.setFieldData(fields_Imag[i],fm);
00153   }
00154 
00155   this->utils.setFieldData(coordVec,fm);
00156   this->utils.setFieldData(weights,fm);
00157   opRegion->postRegistrationSetup(fm);
00158   PHAL::SeparableScatterScalarResponse<EvalT,Traits>::postRegistrationSetup(d,fm);
00159 }
00160 
00161 // **********************************************************************
00162 template<typename EvalT, typename Traits>
00163 void QCAD::ResponseFieldIntegral<EvalT, Traits>::
00164 preEvaluate(typename Traits::PreEvalData workset)
00165 {
00166   for (typename PHX::MDField<ScalarT>::size_type i=0; 
00167        i<this->global_response.size(); i++)
00168     this->global_response[i] = 0.0;
00169 
00170   // Do global initialization
00171   PHAL::SeparableScatterScalarResponse<EvalT,Traits>::preEvaluate(workset);
00172 }
00173 
00174 // **********************************************************************
00175 template<typename EvalT, typename Traits>
00176 void QCAD::ResponseFieldIntegral<EvalT, Traits>::
00177 evaluateFields(typename Traits::EvalData workset)
00178 {
00179   // Zero out local response
00180   for (typename PHX::MDField<ScalarT>::size_type i=0; 
00181        i<this->local_response.size(); i++)
00182     this->local_response[i] = 0.0;
00183 
00184   typename std::vector<PHX::MDField<ScalarT,Cell,QuadPoint> >::const_iterator it;
00185 
00186   if(opRegion->elementBlockIsInRegion(workset.EBName)) {
00187 
00188     ScalarT term, val; //, dbI = 0.0;
00189     std::size_t n, max, nExtraMinuses, nOneBits, nBits = fields.size();
00190 
00191     //DEBUG
00192     //std::size_t nContrib1 = 0, nContrib2 = 0;
00193     //ScalarT dbMaxRe[10], dbMaxIm[10];
00194     //for(std::size_t i=0; i<10; i++) dbMaxRe[i] = dbMaxIm[i] = 0.0;
00195 
00196     for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00197       if(!opRegion->cellIsInRegion(cell)) continue;
00198 
00199       for (std::size_t qp=0; qp < numQPs; ++qp) {
00200   val = 0.0;
00201 
00202   //Loop over all possible combinations of Re/Im parts which form product terms and 
00203   // add the relevant ones (depending on whether we're returning the overall real or
00204   // imaginary part of the integral) to get the integrand value for this (cell,qp).
00205   // We do this by mapping the Re/Im choice onto a string of N bits, where N is the 
00206   // number of fields being multiplied together. (0 = RePart, 1 = ImPart)
00207 
00208   //nContrib1++; //DEBUG
00209 
00210   //for(it = fields.begin(); it != fields.end(); ++it)
00211   max = (std::size_t)std::pow(2.,(int)nBits);
00212 //  max = pow(2.0,static_cast<int>(nBits));
00213   for(std::size_t i=0; i<max; i++) {
00214 
00215     // Count the number of 1 bits, and exit early if 
00216     //  there's a 1 bit for a field that is not complex
00217     nOneBits = nExtraMinuses = 0;
00218     for(n=0; n<nBits; n++) {
00219       if( (0x1 << n) & i ) { // if n-th bit of i is set (use Im part of n-th field)
00220         if(!fieldIsComplex[n]) break;
00221         if(conjugateFieldFlag[n]) nExtraMinuses++;
00222         nOneBits++;
00223       } 
00224     }
00225     if(n < nBits) continue;  // we exited early, signaling this product can't contribute
00226 
00227     //check if this combination of Re/Im parts contributes to the overall Re or Im part we return
00228     if( (bReturnImagPart && nOneBits % 2) || (!bReturnImagPart && nOneBits % 2 == 0)) {
00229       term = (nOneBits % 4 >= 2) ? -1.0 : 1.0; //apply minus sign if nOneBits % 4 == 2 (-1) or == 3 (-i)
00230       if(nExtraMinuses % 2) term *= -1.0;      //apply minus sign due to conjugations
00231       //nContrib2++;
00232 
00233       //multiply fields together
00234       for(std::size_t m=0; m<nBits; m++) { 
00235         if( (0x1 << m) & i ) {
00236     term *= fields_Imag[m](cell,qp);
00237     //if( abs(fields_Imag[m](cell,qp)) > dbMaxIm[m]) dbMaxIm[m] = abs(fields_Imag[m](cell,qp));
00238         } 
00239         else {
00240     term *= fields[m](cell,qp);
00241     //if( abs(fields[m](cell,qp)) > dbMaxRe[m]) dbMaxRe[m] = abs(fields[m](cell,qp));
00242         }
00243       }
00244 
00245       val += term;  //add term to overall integrand
00246     }
00247   }
00248   val *= weights(cell,qp) * scaling; //multiply integrand by volume
00249 
00250   //dbI += val; //DEBUG
00251         this->local_response(cell) += val;
00252   this->global_response(0) += val;
00253       }
00254     }
00255 
00256     //DEBUG
00257     /*if(fieldNames.size() > 1) {
00258       std::cout << "DB: " << (bReturnImagPart == true ? "Im" : "Re") << " Field Integral - int(";
00259       for(std::size_t i=0; i<fieldNames.size(); i++) 
00260   std::cout << fieldNames[i] << "," << (conjugateFieldFlag[i] ? "-" : "") << (fieldIsComplex[i] ? fieldNames_Imag[i] : "X") << " * ";
00261       std::cout << " dV) -- I += " << dbI << "  (ebName = " << workset.EBName << 
00262   " contrib1=" << nContrib1 << " contrib2=" << nContrib2 << ")" << std::endl;
00263       std::cout << "DB MAX of Fields Re: " << dbMaxRe[0] << "," << dbMaxRe[1] << "," << dbMaxRe[2] << "," << dbMaxRe[3] << std::endl;
00264       std::cout << "DB MAX of Fields Im: " << dbMaxIm[0] << "," << dbMaxIm[1] << "," << dbMaxIm[2] << "," << dbMaxIm[3] << std::endl;
00265       }*/
00266   }
00267 
00268   // Do any local-scattering necessary
00269   PHAL::SeparableScatterScalarResponse<EvalT,Traits>::evaluateFields(workset);
00270 }
00271 
00272 // **********************************************************************
00273 template<typename EvalT, typename Traits>
00274 void QCAD::ResponseFieldIntegral<EvalT, Traits>::
00275 postEvaluate(typename Traits::PostEvalData workset)
00276 {
00277   // Add contributions across processors
00278   Teuchos::RCP< Teuchos::ValueTypeSerializer<int,ScalarT> > serializer =
00279     workset.serializerManager.template getValue<EvalT>();
00280   Teuchos::reduceAll(
00281     *workset.comm, *serializer, Teuchos::REDUCE_SUM,
00282     this->global_response.size(), &this->global_response[0], 
00283     &this->global_response[0]);
00284 
00285   if (bPositiveOnly && this->global_response[0] < 1e-6) {
00286     this->global_response[0] = 1e+100;
00287   }
00288   
00289   // Do global scattering
00290   PHAL::SeparableScatterScalarResponse<EvalT,Traits>::postEvaluate(workset);
00291 }
00292 
00293 // **********************************************************************
00294 template<typename EvalT,typename Traits>
00295 Teuchos::RCP<const Teuchos::ParameterList>
00296 QCAD::ResponseFieldIntegral<EvalT,Traits>::getValidResponseParameters() const
00297 {
00298   Teuchos::RCP<Teuchos::ParameterList> validPL =
00299       rcp(new Teuchos::ParameterList("Valid ResponseFieldIntegral Params"));;
00300   Teuchos::RCP<const Teuchos::ParameterList> baseValidPL =
00301     PHAL::SeparableScatterScalarResponse<EvalT,Traits>::getValidResponseParameters();
00302   validPL->setParameters(*baseValidPL);
00303 
00304   Teuchos::RCP<const Teuchos::ParameterList> regionValidPL =
00305     QCAD::MeshRegion<EvalT,Traits>::getValidParameters();
00306   validPL->setParameters(*regionValidPL);
00307 
00308   validPL->set<std::string>("Name", "", "Name of response function");
00309   validPL->set<int>("Phalanx Graph Visualization Detail", 0, "Make dot file to visualize phalanx graph");
00310 
00311   validPL->set<std::string>("Field Name", "", "Name of Field to integrate");
00312   validPL->set<std::string>("Field Name Im", "", "Name of Field to integrate");
00313   validPL->set<bool>("Conjugate Field", false, "Whether a (complex-valued) field should be conjugated in product of fields");
00314   for(int i=1; i < QCAD::MAX_FIELDNAMES_IN_INTEGRAL; i++)
00315     validPL->set<std::string>(Albany::strint("Field Name",i), "", "Name of Field to integrate (multiplied into integrand)");
00316   for(int i=1; i < QCAD::MAX_FIELDNAMES_IN_INTEGRAL; i++)
00317     validPL->set<std::string>(Albany::strint("Field Name Im",i), "", "Name of Imaginar part of Field to integrate (multiplied into integrand)");
00318   for(int i=1; i < QCAD::MAX_FIELDNAMES_IN_INTEGRAL; i++)
00319     validPL->set<bool>(Albany::strint("Conjugate Field",i), false, "Whether field should be conjugated in product of fields");
00320 
00321   validPL->set<std::string>("Integrand Length Unit","cm","Linear length unit of integrand, e.g. cm for integrand in cm^-3.  Can be m, cm, um, nm, or mesh.");
00322   validPL->set<bool>("Positive Return Only",false);
00323   validPL->set<bool>("Return Imaginary Part",false,"True return imaginary part of integral, False returns real part");
00324 
00325 
00326   validPL->set<std::string>("Description", "", "Description of this response used by post processors");
00327   
00328   return validPL;
00329 }
00330 
00331 // **********************************************************************
00332
QCAD_ResponseFieldIntegral_Def.hpp
