Albany: a Trilinos-based PDE code: HDiffusionDeformationMatterResidual

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 <Teuchos_TestForException.hpp>
00008 #include <Phalanx_DataLayout.hpp>
00009 
00010 #include <Intrepid_FunctionSpaceTools.hpp>
00011 //#include <Intrepid_RealSpaceTools.hpp>
00012 #include <Intrepid_MiniTensor.h>
00013 
00014 #include <typeinfo>
00015 
00016 namespace LCM {
00017 
00018   //**********************************************************************
00019   template<typename EvalT, typename Traits>
00020   HDiffusionDeformationMatterResidual<EvalT, Traits>::
00021   HDiffusionDeformationMatterResidual(const Teuchos::ParameterList& p) :
00022     wBF         (p.get<std::string>           ("Weighted BF Name"),
00023                  p.get<Teuchos::RCP<PHX::DataLayout> >("Node QP Scalar Data Layout") ),
00024     wGradBF     (p.get<std::string>           ("Weighted Gradient BF Name"),
00025                  p.get<Teuchos::RCP<PHX::DataLayout> >("Node QP Vector Data Layout") ),
00026     GradBF      (p.get<std::string>           ("Gradient BF Name"),
00027                  p.get<Teuchos::RCP<PHX::DataLayout> >("Node QP Vector Data Layout") ),
00028     Dstar (p.get<std::string>                 ("Effective Diffusivity Name"),
00029            p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00030     DL   (p.get<std::string>                  ("Diffusion Coefficient Name"),
00031           p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00032     Clattice (p.get<std::string>              ("QP Variable Name"),
00033               p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00034     eqps (p.get<std::string>                  ("eqps Name"),
00035           p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00036     eqpsFactor (p.get<std::string>            ("Strain Rate Factor Name"),
00037                 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00038     Ctrapped (p.get<std::string>              ("Trapped Concentration Name"),
00039               p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00040     Ntrapped (p.get<std::string>              ("Trapped Solvent Name"),
00041               p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00042     CLGrad       (p.get<std::string>          ("Gradient QP Variable Name"),
00043                   p.get<Teuchos::RCP<PHX::DataLayout> >("QP Vector Data Layout") ),
00044     stressGrad       (p.get<std::string>      ("Gradient Hydrostatic Stress Name"),
00045                       p.get<Teuchos::RCP<PHX::DataLayout> >("QP Vector Data Layout") ),
00046     DefGrad      (p.get<std::string>          ("Deformation Gradient Name"),
00047                   p.get<Teuchos::RCP<PHX::DataLayout> >("QP Tensor Data Layout") ),
00048     Pstress      (p.get<std::string>          ("Stress Name"),
00049                   p.get<Teuchos::RCP<PHX::DataLayout> >("QP Tensor Data Layout") ),
00050     weights       (p.get<std::string>         ("Weights Name"),
00051                    p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00052     tauFactor  (p.get<std::string>            ("Tau Contribution Name"),
00053                 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00054     elementLength (p.get<std::string>         ("Element Length Name"),
00055                    p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00056     deltaTime (p.get<std::string>             ("Delta Time Name"),
00057                p.get<Teuchos::RCP<PHX::DataLayout> >("Workset Scalar Data Layout")),
00058     TResidual   (p.get<std::string>           ("Residual Name"),
00059                  p.get<Teuchos::RCP<PHX::DataLayout> >("Node Scalar Data Layout") ),
00060      stab_param_(p.get<RealType>("Stabilization Parameter"))
00061 
00062   {
00063     if (p.isType<bool>("Disable Transient"))
00064       enableTransient = !p.get<bool>("Disable Transient");
00065     else enableTransient = true;
00066 
00067     this->addDependentField(elementLength);
00068     this->addDependentField(wBF);
00069     this->addDependentField(wGradBF);
00070     this->addDependentField(GradBF);
00071     this->addDependentField(Dstar);
00072     this->addDependentField(DL);
00073     this->addDependentField(Clattice);
00074     this->addDependentField(eqps);
00075     this->addDependentField(eqpsFactor);
00076     this->addDependentField(Ctrapped);
00077     this->addDependentField(Ntrapped);
00078     this->addDependentField(CLGrad);
00079     this->addDependentField(stressGrad);
00080     this->addDependentField(DefGrad);
00081     this->addDependentField(Pstress);
00082     this->addDependentField(weights);
00083     this->addDependentField(tauFactor);
00084     this->addDependentField(deltaTime);
00085 
00086     this->addEvaluatedField(TResidual);
00087 
00088 
00089     Teuchos::RCP<PHX::DataLayout> vector_dl =
00090       p.get< Teuchos::RCP<PHX::DataLayout> >("QP Vector Data Layout");
00091     std::vector<PHX::DataLayout::size_type> dims;
00092     vector_dl->dimensions(dims);
00093     numQPs  = dims[1];
00094     numDims = dims[2];
00095 
00096     Teuchos::RCP<PHX::DataLayout> node_dl =
00097       p.get< Teuchos::RCP<PHX::DataLayout> >("Node Scalar Data Layout");
00098     std::vector<PHX::DataLayout::size_type> ndims;
00099     node_dl->dimensions(ndims);
00100     worksetSize = dims[0];
00101     numNodes = dims[1];
00102 
00103     // Get data from previous converged time step
00104     ClatticeName = p.get<std::string>("QP Variable Name")+"_old";
00105     CLGradName = p.get<std::string>("Gradient QP Variable Name")+"_old";
00106     eqpsName = p.get<std::string>("eqps Name")+"_old";
00107 
00108     // Allocate workspace for temporary variables
00109     Hflux.resize(worksetSize, numQPs, numDims);
00110     pterm.resize(worksetSize, numQPs);
00111     tpterm.resize(worksetSize, numNodes, numQPs);
00112     artificalDL.resize(worksetSize, numQPs);
00113     stabilizedDL.resize(worksetSize, numQPs);
00114     C.resize(worksetSize, numQPs, numDims, numDims);
00115     Cinv.resize(worksetSize, numQPs, numDims, numDims);
00116     CinvTgrad.resize(worksetSize, numQPs, numDims);
00117     CinvTgrad_old.resize(worksetSize, numQPs, numDims);
00118     CinvTaugrad.resize(worksetSize, numQPs, numDims);
00119 
00120     this->setName("HDiffusionDeformationMatterResidual"+PHX::TypeString<EvalT>::value);
00121 
00122   }
00123 
00124   //**********************************************************************
00125   template<typename EvalT, typename Traits>
00126   void HDiffusionDeformationMatterResidual<EvalT, Traits>::
00127   postRegistrationSetup(typename Traits::SetupData d,
00128                         PHX::FieldManager<Traits>& fm)
00129   {
00130     this->utils.setFieldData(elementLength,fm);
00131     this->utils.setFieldData(wBF,fm);
00132     this->utils.setFieldData(wGradBF,fm);
00133     this->utils.setFieldData(GradBF,fm);
00134     this->utils.setFieldData(Dstar,fm);
00135     this->utils.setFieldData(DL,fm);
00136     this->utils.setFieldData(Clattice,fm);
00137     this->utils.setFieldData(eqps,fm);
00138     this->utils.setFieldData(eqpsFactor,fm);
00139     this->utils.setFieldData(Ctrapped,fm);
00140     this->utils.setFieldData(Ntrapped,fm);
00141     this->utils.setFieldData(CLGrad,fm);
00142     this->utils.setFieldData(stressGrad,fm);
00143     this->utils.setFieldData(DefGrad,fm);
00144     this->utils.setFieldData(Pstress,fm);
00145     this->utils.setFieldData(tauFactor,fm);
00146     this->utils.setFieldData(weights,fm);
00147     this->utils.setFieldData(deltaTime,fm);
00148 
00149     //    if (haveSource) this->utils.setFieldData(Source);
00150     //   if (haveMechSource) this->utils.setFieldData(MechSource);
00151 
00152     this->utils.setFieldData(TResidual,fm);
00153   }
00154 
00155   //**********************************************************************
00156   template<typename EvalT, typename Traits>
00157   void HDiffusionDeformationMatterResidual<EvalT, Traits>::
00158   evaluateFields(typename Traits::EvalData workset)
00159   {
00160    typedef Intrepid::FunctionSpaceTools FST;
00161 //   typedef Intrepid::RealSpaceTools<ScalarT> RST;
00162 
00163     Albany::MDArray Clattice_old = (*workset.stateArrayPtr)[ClatticeName];
00164     Albany::MDArray eqps_old = (*workset.stateArrayPtr)[eqpsName];
00165     Albany::MDArray CLGrad_old = (*workset.stateArrayPtr)[CLGradName];
00166 
00167     ScalarT dt = deltaTime(0);
00168     ScalarT temp(0.0);
00169 
00170 
00171     // compute artifical diffusivity
00172     // for 1D this is identical to lumped mass as shown in Prevost's paper.
00173     for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00174       for (std::size_t qp=0; qp < numQPs; ++qp) {
00175       if (dt == 0){
00176                artificalDL(cell,qp) = 0;
00177       } else {
00178         temp = elementLength(cell,qp)*elementLength(cell,qp)/6.0*Dstar(cell,qp)/DL(cell,qp)/dt;
00179         artificalDL(cell,qp) = stab_param_*
00180           //  (temp) // temp - DL is closer to the limit ...if lumped mass is preferred..
00181           std::abs(temp) // should be 1 but use 0.5 for safety
00182           *(0.5 + 0.5*std::tanh( (temp-1)/DL(cell,qp)  ))
00183           *DL(cell,qp);
00184       }
00185         stabilizedDL(cell,qp) = artificalDL(cell,qp)/( DL(cell,qp) + artificalDL(cell,qp) );
00186       }
00187     }
00188 
00189     // compute the 'material' flux
00190     for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00191 
00192       for (std::size_t qp=0; qp < numQPs; ++qp) {
00193 
00194 
00195         Intrepid::Tensor<ScalarT> F(numDims, &DefGrad(cell, qp, 0, 0));
00196         Intrepid::Tensor<ScalarT> C_tensor_ = Intrepid::t_dot(F,F);
00197         Intrepid::Tensor<ScalarT> C_inv_tensor_ = Intrepid::inverse(C_tensor_);
00198 
00199           Intrepid::Vector<ScalarT> C_grad_(numDims, &CLGrad(cell, qp, 0));
00200           Intrepid::Vector<ScalarT> C_grad_in_ref_ = Intrepid::dot(C_inv_tensor_, C_grad_ );
00201 
00202          for (std::size_t j=0; j<numDims; j++){
00203            Hflux(cell,qp,j) = (1.0 -stabilizedDL(cell,qp))*C_grad_in_ref_(j)*dt;
00204         }
00205       }
00206     }
00207 
00208     FST::integrate<ScalarT>(TResidual, Hflux, wGradBF, Intrepid::COMP_CPP, false); // this also works
00209 
00210     for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00211       for (std::size_t node=0; node < numNodes; ++node) {
00212         for (std::size_t qp=0; qp < numQPs; ++qp) {
00213 
00214                    // Divide the equation by DL to avoid ill-conditioned tangent
00215                    temp =  1.0/ ( DL(cell,qp)  + artificalDL(cell,qp)  );
00216 
00217                   // Transient Term
00218                   TResidual(cell,node) +=  Dstar(cell, qp)*
00219                                                          (Clattice(cell,qp)- Clattice_old(cell, qp) )*
00220                                                              wBF(cell, node, qp)*temp;
00221 
00222                  // Strain Rate Term
00223                  TResidual(cell,node) +=  eqpsFactor(cell,qp)*
00224                                                            (eqps(cell,qp)- eqps_old(cell, qp))*
00225                                                             wBF(cell, node, qp)*temp;
00226 
00227                  // hydrostatic stress term
00228                  for (std::size_t dim=0; dim < numDims; ++dim) {
00229                          TResidual(cell,node) -= tauFactor(cell,qp)*Clattice(cell,qp)*
00230                                                                   wGradBF(cell, node, qp, dim)*
00231                                                                   stressGrad(cell, qp, dim)*dt*temp;
00232                  }
00233             }
00234          }
00235      }
00236 
00237     //---------------------------------------------------------------------------//
00238     // Stabilization Term
00239     ScalarT CLPbar(0);
00240     ScalarT vol(0);
00241 
00242     for (std::size_t cell=0; cell < workset.numCells; ++cell){
00243       CLPbar = 0.0;
00244       vol = 0.0;
00245       for (std::size_t qp=0; qp < numQPs; ++qp) {
00246         CLPbar += weights(cell,qp)*
00247                            (Clattice(cell,qp) - Clattice_old(cell, qp)  );
00248         vol  += weights(cell,qp);
00249       }
00250       CLPbar /= vol;
00251 
00252       for (std::size_t qp=0; qp < numQPs; ++qp) {
00253         pterm(cell,qp) = CLPbar;
00254       }
00255 
00256       for (std::size_t node=0; node < numNodes; ++node) {
00257         trialPbar = 0.0;
00258         for (std::size_t qp=0; qp < numQPs; ++qp) {
00259           trialPbar += wBF(cell,node,qp);
00260         }
00261         trialPbar /= vol;
00262         for (std::size_t qp=0; qp < numQPs; ++qp) {
00263           tpterm(cell,node,qp) = trialPbar;
00264         }
00265       }
00266     }
00267 
00268     for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00269       for (std::size_t node=0; node < numNodes; ++node) {
00270         for (std::size_t qp=0; qp < numQPs; ++qp) {
00271           temp =  1.0/ ( DL(cell,qp)  + artificalDL(cell,qp)  );
00272           TResidual(cell,node) -=  stab_param_*Dstar(cell, qp)*temp*
00273                                                    (-Clattice(cell,qp) + Clattice_old(cell, qp)+pterm(cell,qp))*
00274                                                     wBF(cell, node, qp);
00275         }
00276       }
00277     }
00278 
00279 
00280 }
00281   //**********************************************************************
00282 }
00283 
00284
HDiffusionDeformationMatterResidual_Def.hpp
