Albany: a Trilinos-based PDE code: QCAD_ResponseFieldValue

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 
00011 // **********************************************************************
00012 // Specialization: Jacobian
00013 // **********************************************************************
00014 
00015 template<typename Traits>
00016 void 
00017 QCAD::FieldValueScatterScalarResponse<PHAL::AlbanyTraits::Jacobian, Traits>::
00018 postEvaluate(typename Traits::PostEvalData workset)
00019 {
00020   // Here we scatter the *global* response
00021   Teuchos::RCP<Epetra_Vector> g = workset.g;
00022   if (g != Teuchos::null)
00023     for (int res = 0; res < this->field_components.size(); res++) {
00024       (*g)[res] = this->global_response[this->field_components[res]].val();
00025   }
00026 
00027   Teuchos::RCP<Epetra_MultiVector> dgdx = workset.dgdx;
00028   Teuchos::RCP<Epetra_MultiVector> dgdxdot = workset.dgdxdot;
00029   Teuchos::RCP<Epetra_MultiVector> overlapped_dgdx = workset.overlapped_dgdx;
00030   Teuchos::RCP<Epetra_MultiVector> overlapped_dgdxdot = 
00031     workset.overlapped_dgdxdot;
00032   Teuchos::RCP<Epetra_MultiVector> dg, overlapped_dg;
00033   if (dgdx != Teuchos::null) {
00034     dg = dgdx;
00035     overlapped_dg = overlapped_dgdx;
00036   }
00037   else {
00038     dg = dgdxdot;
00039     overlapped_dg = overlapped_dgdxdot;
00040   }
00041 
00042   dg->PutScalar(0.0);
00043   overlapped_dg->PutScalar(0.0);
00044    
00045   // Extract derivatives for the cell corresponding to nodeID
00046   if (nodeID != Teuchos::null) {
00047 
00048     // Loop over responses
00049     for (int res = 0; res < this->field_components.size(); res++) {
00050       ScalarT& val = this->global_response(this->field_components[res]);
00051       
00052       // Loop over nodes in cell
00053       for (int node_dof=0; node_dof<numNodes; node_dof++) {
00054   int neq = nodeID[node_dof].size();
00055   
00056   // Loop over equations per node
00057   for (int eq_dof=0; eq_dof<neq; eq_dof++) {
00058     
00059     // local derivative component
00060     int deriv = neq * node_dof + eq_dof;
00061     
00062     // local DOF
00063     int dof = nodeID[node_dof][eq_dof];
00064     
00065     // Set dg/dx
00066     overlapped_dg->ReplaceMyValue(dof, res, val.dx(deriv));
00067     
00068   } // column equations
00069       } // column nodes
00070     } // response
00071   } // cell belongs to this proc
00072 
00073   dg->Export(*overlapped_dg, *workset.x_importer, Add);
00074 }
00075 
00076 // **********************************************************************
00077 // Specialization: Stochastic Galerkin Jacobian
00078 // **********************************************************************
00079 
00080 #ifdef ALBANY_SG_MP
00081 template<typename Traits>
00082 void 
00083 QCAD::FieldValueScatterScalarResponse<PHAL::AlbanyTraits::SGJacobian, Traits>::
00084 postEvaluate(typename Traits::PostEvalData workset)
00085 {
00086   // Here we scatter the *global* SG response
00087   Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > g_sg = workset.sg_g;
00088   if (g_sg != Teuchos::null) {
00089     for (int res = 0; res < this->field_components.size(); res++) {
00090       ScalarT& val = this->global_response[this->field_components[res]];
00091       for (int block=0; block<g_sg->size(); block++)
00092   (*g_sg)[block][res] = val.val().coeff(block);
00093     }
00094   }
00095 
00096   // Here we scatter the *global* SG response derivatives
00097   Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> dgdx_sg = 
00098     workset.sg_dgdx;
00099   Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> dgdxdot_sg = 
00100     workset.sg_dgdxdot;
00101   Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> overlapped_dgdx_sg = 
00102     workset.overlapped_sg_dgdx;
00103   Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> overlapped_dgdxdot_sg = 
00104     workset.overlapped_sg_dgdxdot;
00105 
00106   Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> dg_sg, overlapped_dg_sg;
00107   if (dgdx_sg != Teuchos::null) {
00108     dg_sg = dgdx_sg;
00109     overlapped_dg_sg = overlapped_dgdx_sg;
00110   }
00111   else {
00112     dg_sg = dgdxdot_sg;
00113     overlapped_dg_sg = overlapped_dgdxdot_sg;
00114   }
00115 
00116   dg_sg->init(0.0);
00117   overlapped_dg_sg->init(0.0);
00118 
00119   // Extract derivatives for the cell corresponding to nodeID
00120   if (nodeID != Teuchos::null) {
00121 
00122     // Loop over responses
00123     for (int res = 0; res < this->field_components.size(); res++) {
00124       ScalarT& val = this->global_response(this->field_components[res]);
00125       
00126       // Loop over nodes in cell
00127       for (int node_dof=0; node_dof<numNodes; node_dof++) {
00128   int neq = nodeID[node_dof].size();
00129   
00130   // Loop over equations per node
00131   for (int eq_dof=0; eq_dof<neq; eq_dof++) {
00132     
00133     // local derivative component
00134     int deriv = neq * node_dof + eq_dof;
00135     
00136     // local DOF
00137     int dof = nodeID[node_dof][eq_dof];
00138     
00139     // Set dg/dx
00140     for (int block=0; block<dg_sg->size(); block++)
00141       (*overlapped_dg_sg)[block].ReplaceMyValue(dof, res, 
00142                   val.dx(deriv).coeff(block));
00143     
00144   } // column equations
00145       } // response
00146     } // node
00147   } // cell belongs to this proc
00148 
00149   for (int block=0; block<dgdx_sg->size(); block++)
00150     (*dg_sg)[block].Export((*overlapped_dg_sg)[block], 
00151          *workset.x_importer, Add);
00152 }
00153 
00154 // **********************************************************************
00155 // Specialization: Multi-point Jacobian
00156 // **********************************************************************
00157 
00158 template<typename Traits>
00159 void 
00160 QCAD::FieldValueScatterScalarResponse<PHAL::AlbanyTraits::MPJacobian, Traits>::
00161 postEvaluate(typename Traits::PostEvalData workset)
00162 {
00163   // Here we scatter the *global* MP response
00164   Teuchos::RCP<Stokhos::ProductEpetraVector> g_mp = workset.mp_g;
00165   if (g_mp != Teuchos::null) {
00166     for (int res = 0; res < this->field_components.size(); res++) {
00167       ScalarT& val = this->global_response[this->field_components[res]];
00168       for (int block=0; block<g_mp->size(); block++)
00169   (*g_mp)[block][res] = val.val().coeff(block);
00170     }
00171   }
00172 
00173   // Here we scatter the *global* MP response derivatives
00174   Teuchos::RCP<Stokhos::ProductEpetraMultiVector> dgdx_mp = 
00175     workset.mp_dgdx;
00176   Teuchos::RCP<Stokhos::ProductEpetraMultiVector> dgdxdot_mp = 
00177     workset.mp_dgdxdot;
00178   Teuchos::RCP<Stokhos::ProductEpetraMultiVector> overlapped_dgdx_mp = 
00179     workset.overlapped_mp_dgdx;
00180   Teuchos::RCP<Stokhos::ProductEpetraMultiVector> overlapped_dgdxdot_mp = 
00181     workset.overlapped_mp_dgdxdot;
00182   Teuchos::RCP<Stokhos::ProductEpetraMultiVector> dg_mp, overlapped_dg_mp;
00183   if (dgdx_mp != Teuchos::null) {
00184     dg_mp = dgdx_mp;
00185     overlapped_dg_mp = overlapped_dgdx_mp;
00186   }
00187   else {
00188     dg_mp = dgdxdot_mp;
00189     overlapped_dg_mp = overlapped_dgdxdot_mp;
00190   }
00191 
00192   dg_mp->init(0.0);
00193   overlapped_dg_mp->init(0.0);
00194 
00195   // Extract derivatives for the cell corresponding to nodeID
00196   if (nodeID != Teuchos::null) {
00197 
00198     // Loop over responses
00199     for (int res = 0; res < this->field_components.size(); res++) {
00200       ScalarT& val = this->global_response(this->field_components[res]);
00201       
00202       // Loop over nodes in cell
00203       for (int node_dof=0; node_dof<numNodes; node_dof++) {
00204   int neq = nodeID[node_dof].size();
00205   
00206   // Loop over equations per node
00207   for (int eq_dof=0; eq_dof<neq; eq_dof++) {
00208     
00209     // local derivative component
00210     int deriv = neq * node_dof + eq_dof;
00211     
00212     // local DOF
00213     int dof = nodeID[node_dof][eq_dof];
00214     
00215     // Set dg/dx
00216     for (int block=0; block<dg_mp->size(); block++)
00217       (*overlapped_dg_mp)[block].ReplaceMyValue(dof, res, 
00218                   val.dx(deriv).coeff(block));
00219     
00220   } // column equations
00221       } // response
00222     } // node
00223   } // cell belongs to this proc
00224 
00225   for (int block=0; block<dgdx_mp->size(); block++)
00226     (*dg_mp)[block].Export((*overlapped_dg_mp)[block], 
00227          *workset.x_importer, Add);
00228 }
00229 #endif //ALBANY_SG_MP
00230 
00231 template<typename EvalT, typename Traits>
00232 QCAD::ResponseFieldValue<EvalT, Traits>::
00233 ResponseFieldValue(Teuchos::ParameterList& p,
00234        const Teuchos::RCP<Albany::Layouts>& dl) :
00235   coordVec("Coord Vec", dl->qp_vector),
00236   weights("Weights", dl->qp_scalar)
00237 {
00238   // get and validate Response parameter list
00239   Teuchos::ParameterList* plist = 
00240     p.get<Teuchos::ParameterList*>("Parameter List");
00241   Teuchos::RCP<const Teuchos::ParameterList> reflist = 
00242     this->getValidResponseParameters();
00243   plist->validateParameters(*reflist,0);
00244 
00246   Teuchos::RCP<Teuchos::ParameterList> paramsFromProblem = 
00247     p.get< Teuchos::RCP<Teuchos::ParameterList> >("Parameters From Problem");
00248 
00249     // Material database (if given)
00250   if(paramsFromProblem != Teuchos::null)
00251     materialDB = paramsFromProblem->get< Teuchos::RCP<QCAD::MaterialDatabase> >("MaterialDB");
00252   else materialDB = Teuchos::null;
00253 
00254   // number of quad points per cell and dimension of space
00255   Teuchos::RCP<PHX::DataLayout> scalar_dl = dl->qp_scalar;
00256   Teuchos::RCP<PHX::DataLayout> vector_dl = dl->qp_vector;
00257   
00258   std::vector<PHX::DataLayout::size_type> dims;
00259   vector_dl->dimensions(dims);
00260   numQPs  = dims[1];
00261   numDims = dims[2];
00262 
00263   opRegion  = Teuchos::rcp( new QCAD::MeshRegion<EvalT, Traits>("Coord Vec","Weights",*plist,materialDB,dl) );
00264 
00265   // User-specified parameters
00266   operation    = plist->get<std::string>("Operation");
00267 
00268   bOpFieldIsVector = false;
00269   if(plist->isParameter("Operation Vector Field Name")) {
00270     opFieldName  = plist->get<std::string>("Operation Vector Field Name");
00271     bOpFieldIsVector = true;
00272   }
00273   else opFieldName  = plist->get<std::string>("Operation Field Name");
00274  
00275   bRetFieldIsVector = false;
00276   if(plist->isParameter("Return Vector Field Name")) {
00277     retFieldName  = plist->get<std::string>("Return Vector Field Name");
00278     bRetFieldIsVector = true;
00279   }
00280   else retFieldName = plist->get<std::string>("Return Field Name", opFieldName);
00281   bReturnOpField = (opFieldName == retFieldName);
00282 
00283   opX = plist->get<bool>("Operate on x-component", true) && (numDims > 0);
00284   opY = plist->get<bool>("Operate on y-component", true) && (numDims > 1);
00285   opZ = plist->get<bool>("Operate on z-component", true) && (numDims > 2);
00286 
00287   // setup operation field and return field (if it's a different field)
00288   if(bOpFieldIsVector) {
00289     PHX::MDField<ScalarT> f(opFieldName, vector_dl); opField = f; }
00290   else {
00291     PHX::MDField<ScalarT> f(opFieldName, scalar_dl); opField = f; }
00292 
00293   if(!bReturnOpField) {
00294     if(bRetFieldIsVector) {
00295       PHX::MDField<ScalarT> f(retFieldName, vector_dl); retField = f; }
00296     else {
00297       PHX::MDField<ScalarT> f(retFieldName, scalar_dl); retField = f; }
00298   }
00299 
00300   // add dependent fields
00301   this->addDependentField(opField);
00302   this->addDependentField(coordVec);
00303   this->addDependentField(weights);
00304   opRegion->addDependentFields(this);
00305   if(!bReturnOpField) this->addDependentField(retField); //when return field is *different* from op field
00306 
00307   // Set sentinal values for max/min problems 
00308   initVals = Teuchos::Array<double>(5, 0.0);
00309   if( operation == "Maximize" ) initVals[1] = -1e200;
00310   else if( operation == "Minimize" ) initVals[1] = 1e100;
00311   else TEUCHOS_TEST_FOR_EXCEPTION (
00312     true, Teuchos::Exceptions::InvalidParameter, std::endl 
00313     << "Error!  Invalid operation type " << operation << std::endl); 
00314 
00315   this->setName(opFieldName+" Response Field Value"+PHX::TypeString<EvalT>::value);
00316 
00317   // Setup scatter evaluator
00318   std::string global_response_name = 
00319     opFieldName + " Global Response Field Value";
00320   //int worksetSize = scalar_dl->dimension(0);
00321   int responseSize = 5;
00322   Teuchos::RCP<PHX::DataLayout> global_response_layout =
00323     Teuchos::rcp(new PHX::MDALayout<Dim>(responseSize));
00324   PHX::Tag<ScalarT> global_response_tag(global_response_name, 
00325           global_response_layout);
00326   p.set("Stand-alone Evaluator", false);
00327   p.set("Global Response Field Tag", global_response_tag);
00328   this->setup(p,dl);
00329 
00330   // Specify which components of response (in this case 0th and 1st) to 
00331   //  scatter derivatives for.
00332   FieldValueScatterScalarResponse<EvalT,Traits>::field_components.resize(2);
00333   FieldValueScatterScalarResponse<EvalT,Traits>::field_components[0] = 0;
00334   FieldValueScatterScalarResponse<EvalT,Traits>::field_components[1] = 1;
00335 }
00336 
00337 // **********************************************************************
00338 template<typename EvalT, typename Traits>
00339 void QCAD::ResponseFieldValue<EvalT, Traits>::
00340 postRegistrationSetup(typename Traits::SetupData d,
00341                       PHX::FieldManager<Traits>& fm)
00342 {
00343   this->utils.setFieldData(opField,fm);
00344   this->utils.setFieldData(coordVec,fm);
00345   this->utils.setFieldData(weights,fm);
00346   if(!bReturnOpField) this->utils.setFieldData(retField,fm);
00347   opRegion->postRegistrationSetup(fm);
00348   QCAD::FieldValueScatterScalarResponse<EvalT,Traits>::postRegistrationSetup(d,fm);
00349 }
00350 
00351 // **********************************************************************
00352 template<typename EvalT, typename Traits>
00353 void QCAD::ResponseFieldValue<EvalT, Traits>::
00354 preEvaluate(typename Traits::PreEvalData workset)
00355 {
00356   for (typename PHX::MDField<ScalarT>::size_type i=0; 
00357        i<this->global_response.size(); i++)
00358     this->global_response[i] = initVals[i];
00359 
00360   // Do global initialization
00361   QCAD::FieldValueScatterScalarResponse<EvalT,Traits>::preEvaluate(workset);
00362 }
00363 
00364 // **********************************************************************
00365 template<typename EvalT, typename Traits>
00366 void QCAD::ResponseFieldValue<EvalT, Traits>::
00367 evaluateFields(typename Traits::EvalData workset)
00368 {
00369   ScalarT opVal, qpVal, cellVol;
00370 
00371   if(!opRegion->elementBlockIsInRegion(workset.EBName))
00372     return;
00373 
00374   for (std::size_t cell=0; cell < workset.numCells; ++cell) 
00375   {
00376     if(!opRegion->cellIsInRegion(cell)) continue;
00377 
00378     // Get the cell volume, used for averaging over a cell
00379     cellVol = 0.0;
00380     for (std::size_t qp=0; qp < numQPs; ++qp)
00381       cellVol += weights(cell,qp);
00382 
00383     // Get the scalar value of the field being operated on which will be used
00384     //  in the operation (all operations just deal with scalar data so far)
00385     opVal = 0.0;
00386     for (std::size_t qp=0; qp < numQPs; ++qp) {
00387       qpVal = 0.0;
00388       if(bOpFieldIsVector) {
00389         if(opX) qpVal += opField(cell,qp,0) * opField(cell,qp,0);
00390         if(opY) qpVal += opField(cell,qp,1) * opField(cell,qp,1);
00391         if(opZ) qpVal += opField(cell,qp,2) * opField(cell,qp,2);
00392       }
00393       else qpVal = opField(cell,qp);
00394       opVal += qpVal * weights(cell,qp);
00395     }
00396     opVal /= cellVol;  
00397     // opVal = the average value of the field operated on over the current cell
00398 
00399       
00400     // Check if the currently stored min/max value needs to be updated
00401     if( (operation == "Maximize" && opVal > this->global_response[1]) ||
00402         (operation == "Minimize" && opVal < this->global_response[1]) ) {
00403       max_nodeID = workset.wsElNodeEqID[cell];
00404 
00405       // set g[0] = value of return field at the current cell (avg)
00406       this->global_response[0]=0.0;
00407       if(bReturnOpField) {
00408         for (std::size_t qp=0; qp < numQPs; ++qp) {
00409           qpVal = 0.0;
00410           if(bOpFieldIsVector) {
00411             for(std::size_t i=0; i<numDims; i++) {
00412               qpVal += opField(cell,qp,i)*opField(cell,qp,i);
00413             }
00414           }
00415           else qpVal = opField(cell,qp);
00416           this->global_response[0] += qpVal * weights(cell,qp);
00417         }
00418       }
00419       else {
00420         for (std::size_t qp=0; qp < numQPs; ++qp) {
00421           qpVal = 0.0;
00422           if(bRetFieldIsVector) {
00423             for(std::size_t i=0; i<numDims; i++) {
00424               qpVal += retField(cell,qp,i)*retField(cell,qp,i);
00425             }
00426           }
00427           else qpVal = retField(cell,qp);
00428           this->global_response[0] += qpVal * weights(cell,qp);
00429         }
00430       }
00431       this->global_response[0] /= cellVol;
00432 
00433       // set g[1] = value of the field operated on at the current cell (avg)
00434       this->global_response[1] = opVal;
00435 
00436       // set g[2+] = average qp coordinate values of the current cell
00437       for(std::size_t i=0; i<numDims; i++) {
00438         this->global_response[i+2] = 0.0;
00439         for (std::size_t qp=0; qp < numQPs; ++qp) 
00440     this->global_response[i+2] += coordVec(cell,qp,i);
00441   this->global_response[i+2] /= numQPs;
00442       }
00443     }
00444 
00445   } // end of loop over cells
00446 
00447   // No local scattering
00448 }
00449 
00450 // **********************************************************************
00451 template<typename EvalT, typename Traits>
00452 void QCAD::ResponseFieldValue<EvalT, Traits>::
00453 postEvaluate(typename Traits::PostEvalData workset)
00454 {
00455   int indexToMax = 1;
00456   ScalarT max = this->global_response[indexToMax];
00457   Teuchos::EReductionType reductType;
00458   if (operation == "Maximize")
00459     reductType = Teuchos::REDUCE_MAX;
00460   else
00461     reductType = Teuchos::REDUCE_MIN;
00462 
00463   Teuchos::RCP< Teuchos::ValueTypeSerializer<int,ScalarT> > serializer =
00464     workset.serializerManager.template getValue<EvalT>();
00465 
00466   // Compute contributions across processors
00467   Teuchos::reduceAll(
00468     *workset.comm, *serializer, reductType, 1, 
00469     &this->global_response[indexToMax], &max);
00470 
00471   int procToBcast;
00472   if( this->global_response[indexToMax] == max ) 
00473     procToBcast = workset.comm->getRank();
00474   else procToBcast = -1;
00475 
00476   int winner;
00477   Teuchos::reduceAll(
00478     *workset.comm, Teuchos::REDUCE_MAX, 1, &procToBcast, &winner);
00479   Teuchos::broadcast(
00480     *workset.comm, *serializer, winner, this->global_response.size(), 
00481     &this->global_response[0]);
00482 
00483   // Do global scattering
00484   if (workset.comm->getRank() == winner)
00485     QCAD::FieldValueScatterScalarResponse<EvalT,Traits>::setNodeID(max_nodeID);
00486   else
00487     QCAD::FieldValueScatterScalarResponse<EvalT,Traits>::setNodeID(Teuchos::null);
00488 
00489   QCAD::FieldValueScatterScalarResponse<EvalT,Traits>::postEvaluate(workset);
00490 }
00491 
00492 // **********************************************************************
00493 template<typename EvalT,typename Traits>
00494 Teuchos::RCP<const Teuchos::ParameterList>
00495 QCAD::ResponseFieldValue<EvalT,Traits>::getValidResponseParameters() const
00496 {
00497   Teuchos::RCP<Teuchos::ParameterList> validPL =
00498       rcp(new Teuchos::ParameterList("Valid ResponseFieldValue Params"));
00499   Teuchos::RCP<const Teuchos::ParameterList> baseValidPL =
00500     QCAD::FieldValueScatterScalarResponse<EvalT,Traits>::getValidResponseParameters();
00501   validPL->setParameters(*baseValidPL);
00502 
00503   Teuchos::RCP<const Teuchos::ParameterList> regionValidPL =
00504     QCAD::MeshRegion<EvalT,Traits>::getValidParameters();
00505   validPL->setParameters(*regionValidPL);
00506 
00507   validPL->set<std::string>("Name", "", "Name of response function");
00508   validPL->set<int>("Phalanx Graph Visualization Detail", 0, "Make dot file to visualize phalanx graph");
00509 
00510   validPL->set<std::string>("Operation", "Maximize", "Operation to perform");
00511   validPL->set<std::string>("Operation Field Name", "", "Scalar field to perform operation on");
00512   validPL->set<std::string>("Operation Vector Field Name", "", "Vector field to perform operation on");
00513   validPL->set<std::string>("Return Field Name", "<operation field name>",
00514            "Scalar field to return value from");
00515   validPL->set<std::string>("Return Vector Field Name", "<operation vector field name>",
00516            "Vector field to return value from");
00517 
00518   validPL->set<bool>("Operate on x-component", true, 
00519          "Whether to perform operation on x component of vector field");
00520   validPL->set<bool>("Operate on y-component", true, 
00521          "Whether to perform operation on y component of vector field");
00522   validPL->set<bool>("Operate on z-component", true, 
00523          "Whether to perform operation on z component of vector field");
00524 
00525   validPL->set<std::string>("Description", "", "Description of this response used by post processors");
00526 
00527   return validPL;
00528 }
00529 
00530 // **********************************************************************
00531
QCAD_ResponseFieldValue_Def.hpp
