Albany: a Trilinos-based PDE code: PHAL

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 #ifndef PHAL_WORKSET_HPP
00008 #define PHAL_WORKSET_HPP
00009 
00010 #include <list>
00011 
00012 #include "Phalanx_ConfigDefs.hpp" // for std::vector
00013 #include "Albany_DataTypes.hpp" 
00014 #include "Epetra_Vector.h"
00015 #include "Epetra_CrsMatrix.h"
00016 #include "Albany_AbstractDiscretization.hpp"
00017 #include "Albany_StateManager.hpp"
00018 #include <Intrepid_FieldContainer.hpp>
00019 
00020 #include "Stokhos_OrthogPolyExpansion.hpp"
00021 #include "Stokhos_EpetraVectorOrthogPoly.hpp"
00022 #include "Stokhos_EpetraMultiVectorOrthogPoly.hpp"
00023 
00024 #include "PHAL_AlbanyTraits.hpp"
00025 #include "PHAL_TypeKeyMap.hpp"
00026 #include "Teuchos_RCP.hpp"
00027 #include "Teuchos_Comm.hpp"
00028 #include "Epetra_Import.h"
00029 
00030 namespace PHAL {
00031 
00032 struct Workset {
00033 
00034   typedef AlbanyTraits::EvalTypes ET;
00035   
00036   Workset() :
00037     transientTerms(false), accelerationTerms(false), ignore_residual(false) {}
00038 
00039   unsigned int numCells;
00040   unsigned int wsIndex;
00041 
00042   Teuchos::RCP<Stokhos::OrthogPolyExpansion<int,double> > sg_expansion;
00043 
00044   Teuchos::RCP<const Epetra_Vector> x;
00045   Teuchos::RCP<const Epetra_Vector> xdot;
00046   Teuchos::RCP<const Epetra_Vector> xdotdot;
00047   Teuchos::RCP<ParamVec> params;
00048   Teuchos::RCP<const Epetra_MultiVector> Vx;
00049   Teuchos::RCP<const Epetra_MultiVector> Vxdot;
00050   Teuchos::RCP<const Epetra_MultiVector> Vxdotdot;
00051   Teuchos::RCP<const Epetra_MultiVector> Vp;
00052   Teuchos::RCP<const Stokhos::EpetraVectorOrthogPoly > sg_x;
00053 
00054   Teuchos::RCP<const Stokhos::EpetraVectorOrthogPoly > sg_xdot;
00055   Teuchos::RCP<const Stokhos::EpetraVectorOrthogPoly > sg_xdotdot;
00056   Teuchos::RCP<const Stokhos::ProductEpetraVector > mp_x;
00057   Teuchos::RCP<const Stokhos::ProductEpetraVector > mp_xdot;
00058   Teuchos::RCP<const Stokhos::ProductEpetraVector > mp_xdotdot;
00059 
00060   Teuchos::RCP<Epetra_Vector> f;
00061   Teuchos::RCP<Epetra_CrsMatrix> Jac;
00062   Teuchos::RCP<Epetra_MultiVector> JV;
00063   Teuchos::RCP<Epetra_MultiVector> fp;
00064   Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > sg_f;
00065   Teuchos::RCP< Stokhos::VectorOrthogPoly<Epetra_CrsMatrix> > sg_Jac;
00066   Teuchos::RCP< Stokhos::EpetraMultiVectorOrthogPoly > sg_JV;
00067   Teuchos::RCP< Stokhos::EpetraMultiVectorOrthogPoly > sg_fp;
00068   Teuchos::RCP< Stokhos::ProductEpetraVector > mp_f;
00069   Teuchos::RCP< Stokhos::ProductContainer<Epetra_CrsMatrix> > mp_Jac;
00070   Teuchos::RCP< Stokhos::ProductEpetraMultiVector > mp_JV;
00071   Teuchos::RCP< Stokhos::ProductEpetraMultiVector > mp_fp;
00072 
00073   Teuchos::RCP<const Albany::NodeSetList> nodeSets;
00074   Teuchos::RCP<const Albany::NodeSetCoordList> nodeSetCoords;
00075 
00076   Teuchos::RCP<const Albany::SideSetList> sideSets;
00077 
00078   // jacobian and mass matrix coefficients for matrix fill
00079   double j_coeff;
00080   double m_coeff; //d(x_dot)/dx_{new}
00081   double n_coeff; //d(x_dotdot)/dx_{new}
00082 
00083   // Current Time as defined by Rythmos
00084   double current_time;
00085   double previous_time;
00086  
00087   // flag indicating whether to sum tangent derivatives, i.e.,
00088   // compute alpha*df/dxdot*Vxdot + beta*df/dx*Vx + omega*df/dxddotot*Vxdotdot + df/dp*Vp or
00089   // compute alpha*df/dxdot*Vxdot + beta*df/dx*Vx + omega*df/dxdotdot*Vxdotdot and df/dp*Vp separately
00090   int num_cols_x;
00091   int num_cols_p;
00092   int param_offset;
00093 
00094   std::vector<int> *coord_deriv_indices;
00095 
00096   Teuchos::ArrayRCP<Teuchos::ArrayRCP<Teuchos::ArrayRCP<int> > >  wsElNodeEqID;
00097   Teuchos::ArrayRCP<Teuchos::ArrayRCP<int> >  wsElNodeID;
00098   Teuchos::ArrayRCP<Teuchos::ArrayRCP<double*> >  wsCoords;
00099   Teuchos::ArrayRCP<Teuchos::ArrayRCP<double> >  wsSHeight;
00100   Teuchos::ArrayRCP<Teuchos::ArrayRCP<double> >  wsVolume; // DJL DEBUGGING
00101   Teuchos::ArrayRCP<double>  wsTemperature;
00102   Teuchos::ArrayRCP<Teuchos::ArrayRCP<double> >  wsBasalFriction;
00103   Teuchos::ArrayRCP<Teuchos::ArrayRCP<double> >  wsThickness;
00104   Teuchos::ArrayRCP<double>  wsFlowFactor;
00105   Teuchos::ArrayRCP<Teuchos::ArrayRCP<double*> > wsSurfaceVelocity;
00106   Teuchos::ArrayRCP<Teuchos::ArrayRCP<double*> > wsVelocityRMS;
00107   Teuchos::ArrayRCP<Teuchos::ArrayRCP<Teuchos::ArrayRCP<Teuchos::ArrayRCP<double> > > >  ws_coord_derivs;
00108   std::string EBName;
00109 
00110   Albany::StateArray* stateArrayPtr;
00111   Teuchos::RCP<Albany::EigendataStruct> eigenDataPtr;
00112   Teuchos::RCP<Epetra_MultiVector> auxDataPtr;
00113 
00114   bool transientTerms;
00115   bool accelerationTerms;
00116 
00117   // Flag indicating whether to ignore residual calculations in the 
00118   // Jacobian calculation.  This only works for some problems where the 
00119   // the calculation of the Jacobian doesn't require calculation of the
00120   // residual (such as linear problems), but if it does work it can
00121   // significantly reduce Jacobian calculation cost.
00122   bool ignore_residual;
00123 
00124   // Flag indicated whether we are solving the adjoint operator or the
00125   // forward operator.  This is used in the Albany application when
00126   // either the Jacobian or the transpose of the Jacobian is scattered. 
00127   bool is_adjoint;
00128 
00129   // New field manager response stuff
00130   Teuchos::RCP<const Teuchos::Comm<int> > comm;
00131   Teuchos::RCP<const Epetra_Import> x_importer;
00132   Teuchos::RCP<Epetra_Vector> g;
00133   Teuchos::RCP<Epetra_MultiVector> dgdx;
00134   Teuchos::RCP<Epetra_MultiVector> dgdxdot;
00135   Teuchos::RCP<Epetra_MultiVector> dgdxdotdot;
00136   Teuchos::RCP<Epetra_MultiVector> overlapped_dgdx;
00137   Teuchos::RCP<Epetra_MultiVector> overlapped_dgdxdot;
00138   Teuchos::RCP<Epetra_MultiVector> overlapped_dgdxdotdot;
00139   Teuchos::RCP<Epetra_MultiVector> dgdp;
00140   Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > sg_g;
00141   Teuchos::RCP< Stokhos::EpetraMultiVectorOrthogPoly > sg_dgdx;
00142   Teuchos::RCP< Stokhos::EpetraMultiVectorOrthogPoly > sg_dgdxdot;
00143   Teuchos::RCP< Stokhos::EpetraMultiVectorOrthogPoly > sg_dgdxdotdot;
00144   Teuchos::RCP< Stokhos::EpetraMultiVectorOrthogPoly > overlapped_sg_dgdx;
00145   Teuchos::RCP< Stokhos::EpetraMultiVectorOrthogPoly > overlapped_sg_dgdxdot;
00146   Teuchos::RCP< Stokhos::EpetraMultiVectorOrthogPoly > overlapped_sg_dgdxdotdot;
00147   Teuchos::RCP< Stokhos::EpetraMultiVectorOrthogPoly > sg_dgdp;
00148   Teuchos::RCP< Stokhos::ProductEpetraVector > mp_g;
00149   Teuchos::RCP< Stokhos::ProductEpetraMultiVector > mp_dgdx;
00150   Teuchos::RCP< Stokhos::ProductEpetraMultiVector > mp_dgdxdot;
00151   Teuchos::RCP< Stokhos::ProductEpetraMultiVector > mp_dgdxdotdot;
00152   Teuchos::RCP< Stokhos::ProductEpetraMultiVector > overlapped_mp_dgdx;
00153   Teuchos::RCP< Stokhos::ProductEpetraMultiVector > overlapped_mp_dgdxdot;
00154   Teuchos::RCP< Stokhos::ProductEpetraMultiVector > overlapped_mp_dgdxdotdot;
00155   Teuchos::RCP< Stokhos::ProductEpetraMultiVector > mp_dgdp;
00156   
00157   // Meta-function class encoding T<EvalT::ScalarT> given EvalT
00158   // where T is any lambda expression (typically a placeholder expression)
00159   template <typename T>
00160   struct ApplyEvalT {
00161     template <typename EvalT> struct apply {
00162       typedef typename boost::mpl::apply<T, typename EvalT::ScalarT>::type type;
00163     };
00164   };
00165 
00166   // Meta-function class encoding RCP<ValueTypeSerializer<int,T> > for a given
00167   // type T.  This is to eliminate an error when using a placeholder expression
00168   // for the same thing in CreateLambdaKeyMap below
00169   struct ApplyVTS {
00170     template <typename T>
00171     struct apply {
00172       typedef Teuchos::RCP< Teuchos::ValueTypeSerializer<int,T> > type;
00173     };
00174   };
00175 
00176   // mpl::vector mapping evaluation type EvalT to serialization class
00177   // ValueTypeSerializer<int, EvalT::ScalarT>, which is used for MPI
00178   // communication of scalar types.
00179   typedef PHAL::CreateLambdaKeyMap<AlbanyTraits::BEvalTypes,
00180            ApplyEvalT<ApplyVTS> >::type SerializerMap;
00181 
00182   // Container storing serializers for each evaluation type
00183   PHAL::TypeKeyMap<SerializerMap> serializerManager;
00184 
00185   void print(std::ostream &os){
00186 
00187     os << "Printing workset data:" << std::endl;
00188     os << "\tEB name : " << EBName << std::endl;
00189     os << "\tnumCells : " << numCells << std::endl;
00190     os << "\twsElNodeEqID : " << std::endl;
00191     for(int i = 0; i < wsElNodeEqID.size(); i++)
00192       for(int j = 0; j < wsElNodeEqID[i].size(); j++)
00193         for(int k = 0; k < wsElNodeEqID[i][j].size(); k++)
00194           os << "\t\twsElNodeEqID[" << i << "][" << j << "][" << k << "] = " << 
00195             wsElNodeEqID[i][j][k] << std::endl;
00196     os << "\twsCoords : " << std::endl;
00197     for(int i = 0; i < wsCoords.size(); i++)
00198       for(int j = 0; j < wsCoords[i].size(); j++)
00199           os << "\t\tcoord0:" << wsCoords[i][j][0] << "][" << wsCoords[i][j][1] << std::endl;
00200   }
00201   
00202 };
00203 
00204   template <typename EvalT> struct BuildSerializer {
00205     BuildSerializer(Workset& workset) {}
00206   };
00207   template <> struct BuildSerializer<PHAL::AlbanyTraits::Residual> {
00208     BuildSerializer(Workset& workset) {
00209       Teuchos::RCP< Teuchos::ValueTypeSerializer<int,RealType> > serializer =
00210   Teuchos::rcp(new Teuchos::ValueTypeSerializer<int,RealType>);
00211       workset.serializerManager. 
00212   setValue<PHAL::AlbanyTraits::Residual>(serializer);
00213     }
00214   };
00215   template <> struct BuildSerializer<PHAL::AlbanyTraits::Jacobian> {
00216     BuildSerializer(Workset& workset) {
00217       int num_nodes = workset.wsElNodeEqID[0].size();
00218       int num_eqns =  workset.wsElNodeEqID[0][0].size();
00219       int num_dof = num_nodes * num_eqns;
00220       Teuchos::RCP< Teuchos::ValueTypeSerializer<int,RealType> > 
00221   real_serializer =
00222   Teuchos::rcp(new Teuchos::ValueTypeSerializer<int,RealType>);
00223       Teuchos::RCP< Teuchos::ValueTypeSerializer<int,FadType> > serializer =
00224   Teuchos::rcp(new Teuchos::ValueTypeSerializer<int,FadType>(
00225            real_serializer, num_dof));
00226       workset.serializerManager. 
00227   setValue<PHAL::AlbanyTraits::Jacobian>(serializer);
00228     }
00229   };
00230   template <> struct BuildSerializer<PHAL::AlbanyTraits::Tangent> {
00231     BuildSerializer(Workset& workset) {
00232       int num_cols_tot = workset.param_offset + workset.num_cols_p;
00233       Teuchos::RCP< Teuchos::ValueTypeSerializer<int,RealType> > 
00234   real_serializer =
00235   Teuchos::rcp(new Teuchos::ValueTypeSerializer<int,RealType>);
00236       Teuchos::RCP< Teuchos::ValueTypeSerializer<int,TanFadType> > serializer =
00237   Teuchos::rcp(new Teuchos::ValueTypeSerializer<int,TanFadType>(
00238            real_serializer, num_cols_tot));
00239       workset.serializerManager. 
00240   setValue<PHAL::AlbanyTraits::Tangent>(serializer);
00241     }
00242   };
00243 #ifdef ALBANY_SG_MP
00244   template <> struct BuildSerializer<PHAL::AlbanyTraits::SGResidual> {
00245     BuildSerializer(Workset& workset) {
00246       Teuchos::RCP< Teuchos::ValueTypeSerializer<int,RealType> > 
00247   real_serializer =
00248   Teuchos::rcp(new Teuchos::ValueTypeSerializer<int,RealType>);
00249       Teuchos::RCP< Teuchos::ValueTypeSerializer<int,SGType> > serializer =
00250   Teuchos::rcp(new Teuchos::ValueTypeSerializer<int,SGType>(
00251            workset.sg_expansion, real_serializer));
00252       workset.serializerManager. 
00253   setValue<PHAL::AlbanyTraits::SGResidual>(serializer);
00254     }
00255   };
00256   template <> struct BuildSerializer<PHAL::AlbanyTraits::SGJacobian> {
00257     BuildSerializer(Workset& workset) {
00258       int num_nodes = workset.wsElNodeEqID[0].size();
00259       int num_eqns =  workset.wsElNodeEqID[0][0].size();
00260       int num_dof = num_nodes * num_eqns;
00261       Teuchos::RCP< Teuchos::ValueTypeSerializer<int,RealType> > 
00262   real_serializer =
00263   Teuchos::rcp(new Teuchos::ValueTypeSerializer<int,RealType>);
00264       Teuchos::RCP< Teuchos::ValueTypeSerializer<int,SGType> > sg_serializer =
00265   Teuchos::rcp(new Teuchos::ValueTypeSerializer<int,SGType>(
00266            workset.sg_expansion, real_serializer));
00267       Teuchos::RCP< Teuchos::ValueTypeSerializer<int,SGFadType> > serializer =
00268   Teuchos::rcp(new Teuchos::ValueTypeSerializer<int,SGFadType>(
00269            sg_serializer, num_dof));
00270       workset.serializerManager. 
00271   setValue<PHAL::AlbanyTraits::SGJacobian>(serializer);
00272     }
00273   };
00274   template <> struct BuildSerializer<PHAL::AlbanyTraits::SGTangent> {
00275     BuildSerializer(Workset& workset) {
00276       int num_cols_tot = workset.param_offset + workset.num_cols_p;
00277       Teuchos::RCP< Teuchos::ValueTypeSerializer<int,RealType> > 
00278   real_serializer =
00279   Teuchos::rcp(new Teuchos::ValueTypeSerializer<int,RealType>);
00280       Teuchos::RCP< Teuchos::ValueTypeSerializer<int,SGType> > sg_serializer =
00281   Teuchos::rcp(new Teuchos::ValueTypeSerializer<int,SGType>(
00282            workset.sg_expansion, real_serializer));
00283       Teuchos::RCP< Teuchos::ValueTypeSerializer<int,SGFadType> > serializer =
00284   Teuchos::rcp(new Teuchos::ValueTypeSerializer<int,SGFadType>(
00285            sg_serializer, num_cols_tot));
00286       workset.serializerManager. 
00287   setValue<PHAL::AlbanyTraits::SGTangent>(serializer);
00288     }
00289   };
00290   template <> struct BuildSerializer<PHAL::AlbanyTraits::MPResidual> {
00291     BuildSerializer(Workset& workset) {
00292       int nblock = workset.mp_x->size();
00293       Teuchos::RCP< Teuchos::ValueTypeSerializer<int,RealType> > 
00294   real_serializer =
00295   Teuchos::rcp(new Teuchos::ValueTypeSerializer<int,RealType>);
00296       Teuchos::RCP< Teuchos::ValueTypeSerializer<int,MPType> > serializer =
00297   Teuchos::rcp(new Teuchos::ValueTypeSerializer<int,MPType>(
00298            real_serializer, nblock));
00299       workset.serializerManager. 
00300   setValue<PHAL::AlbanyTraits::MPResidual>(serializer);
00301     }
00302   };
00303   template <> struct BuildSerializer<PHAL::AlbanyTraits::MPJacobian> {
00304     BuildSerializer(Workset& workset) {
00305       int nblock = workset.mp_x->size();
00306       int num_nodes = workset.wsElNodeEqID[0].size();
00307       int num_eqns =  workset.wsElNodeEqID[0][0].size();
00308       int num_dof = num_nodes * num_eqns;
00309        Teuchos::RCP< Teuchos::ValueTypeSerializer<int,RealType> > 
00310    real_serializer =
00311    Teuchos::rcp(new Teuchos::ValueTypeSerializer<int,RealType>);
00312        Teuchos::RCP< Teuchos::ValueTypeSerializer<int,MPType> > mp_serializer =
00313    Teuchos::rcp(new Teuchos::ValueTypeSerializer<int,MPType>(
00314       real_serializer, nblock));
00315        Teuchos::RCP< Teuchos::ValueTypeSerializer<int,MPFadType> > serializer =
00316    Teuchos::rcp(new Teuchos::ValueTypeSerializer<int,MPFadType>(
00317       mp_serializer, num_dof));
00318        workset.serializerManager. 
00319    setValue<PHAL::AlbanyTraits::MPJacobian>(serializer);
00320     }
00321   };
00322   template <> struct BuildSerializer<PHAL::AlbanyTraits::MPTangent> {
00323     BuildSerializer(Workset& workset) {
00324       int nblock = workset.mp_x->size();
00325       int num_cols_tot = workset.param_offset + workset.num_cols_p;
00326       Teuchos::RCP< Teuchos::ValueTypeSerializer<int,RealType> > 
00327   real_serializer =
00328   Teuchos::rcp(new Teuchos::ValueTypeSerializer<int,RealType>);
00329       Teuchos::RCP< Teuchos::ValueTypeSerializer<int,MPType> > mp_serializer =
00330   Teuchos::rcp(new Teuchos::ValueTypeSerializer<int,MPType>(
00331            real_serializer, nblock));
00332       Teuchos::RCP< Teuchos::ValueTypeSerializer<int,MPFadType> > serializer =
00333   Teuchos::rcp(new Teuchos::ValueTypeSerializer<int,MPFadType>(
00334            mp_serializer, num_cols_tot));
00335       workset.serializerManager. 
00336   setValue<PHAL::AlbanyTraits::MPTangent>(serializer);
00337     }
00338   };
00339 #endif //ALBANY_SG_MP
00340 
00341 }
00342 
00343 #endif
PHAL_Workset.hpp
