Albany: a Trilinos-based PDE code: QCAD_SaddleValueResponseFunction.hpp 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 #ifndef QCAD_SADDLEVALUERESPONSEFUNCTION_HPP
00008 #define QCAD_SADDLEVALUERESPONSEFUNCTION_HPP
00009 
00010 #include "Albany_FieldManagerScalarResponseFunction.hpp"
00011 #include "QCAD_MaterialDatabase.hpp"
00012 
00013 #define MAX_DIMENSIONS 3
00014 
00015 namespace QCAD {
00016 
00017   // Helper class: a vector with math operators
00018   class mathVector
00019   {
00020   public:
00021     mathVector();
00022     mathVector(int n);
00023     mathVector(const mathVector& copy);
00024     ~mathVector();
00025 
00026     void resize(std::size_t n);
00027     void fill(double d);
00028     void fill(const double* vec);
00029     double dot(const mathVector& v2) const;
00030     double distanceTo(const mathVector& v2) const;
00031     double distanceTo(const double* p) const;
00032 
00033     double norm() const;
00034     double norm2() const;
00035     void normalize();
00036 
00037     double* data();
00038     const double* data() const;
00039     std::size_t size() const;
00040 
00041     mathVector& operator=(const mathVector& rhs);
00042 
00043     mathVector operator+(const mathVector& v2) const;
00044     mathVector operator-(const mathVector& v2) const;
00045     mathVector operator*(double scale) const;
00046 
00047     mathVector& operator+=(const mathVector& v2);
00048     mathVector& operator-=(const mathVector& v2);
00049     mathVector& operator*=(double scale);
00050     mathVector& operator/=(double scale);
00051 
00052     double& operator[](int i);
00053     const double& operator[](int i) const;
00054 
00055   private:
00056     int dim_;
00057     std::vector<double> data_;
00058   };
00059 
00060   // Data Structure for an image point
00061   struct nebImagePt {
00062     void init(int nDims, double ptRadius) {
00063       coords.resize(nDims); coords.fill(0.0);
00064       velocity.resize(nDims); velocity.fill(0.0);
00065       grad.resize(nDims); grad.fill(0.0);
00066       value = weight = 0.0;
00067       radius = ptRadius;
00068     }
00069 
00070     void init(const mathVector& coordPt, double ptRadius) {
00071       init(coordPt.size(), ptRadius);
00072       coords = coordPt;
00073     }
00074       
00075     mathVector coords;
00076     mathVector velocity;
00077     mathVector grad;
00078     double value;
00079     double weight;
00080     double radius;
00081   };
00082 
00083   std::ostream& operator<<(std::ostream& os, const mathVector& mv);
00084   std::ostream& operator<<(std::ostream& os, const nebImagePt& np);
00085 
00086   // a double array with maximal dimension
00087   struct maxDimPt { 
00088     maxDimPt(const double* p, int numDims) {
00089       for(int i=0; i<numDims; i++) data[i] = p[i];
00090     }
00091     double data[MAX_DIMENSIONS];
00092   };
00093 
00094  
00098   class SaddleValueResponseFunction : 
00099     public Albany::FieldManagerScalarResponseFunction {
00100   public:
00101   
00103     SaddleValueResponseFunction(
00104       const Teuchos::RCP<Albany::Application>& application,
00105       const Teuchos::RCP<Albany::AbstractProblem>& problem,
00106       const Teuchos::RCP<Albany::MeshSpecsStruct>&  ms,
00107       const Teuchos::RCP<Albany::StateManager>& stateMgr,
00108       Teuchos::ParameterList& responseParams);
00109 
00111     virtual ~SaddleValueResponseFunction();
00112 
00114     virtual unsigned int numResponses() const;
00115 
00117     virtual void 
00118     evaluateResponse(const double current_time,
00119          const Epetra_Vector* xdot,
00120          const Epetra_Vector* xdotdot,
00121          const Epetra_Vector& x,
00122          const Teuchos::Array<ParamVec>& p,
00123          Epetra_Vector& g);
00124 
00126     virtual void 
00127     evaluateTangent(const double alpha, 
00128         const double beta,
00129         const double omega,
00130         const double current_time,
00131         bool sum_derivs,
00132         const Epetra_Vector* xdot,
00133         const Epetra_Vector* xdotdot,
00134         const Epetra_Vector& x,
00135         const Teuchos::Array<ParamVec>& p,
00136         ParamVec* deriv_p,
00137         const Epetra_MultiVector* Vxdot,
00138         const Epetra_MultiVector* Vxdotdot,
00139         const Epetra_MultiVector* Vx,
00140         const Epetra_MultiVector* Vp,
00141         Epetra_Vector* g,
00142         Epetra_MultiVector* gx,
00143         Epetra_MultiVector* gp);
00144 
00146     virtual void 
00147     evaluateGradient(const double current_time,
00148          const Epetra_Vector* xdot,
00149          const Epetra_Vector* xdotdot,
00150          const Epetra_Vector& x,
00151          const Teuchos::Array<ParamVec>& p,
00152          ParamVec* deriv_p,
00153          Epetra_Vector* g,
00154          Epetra_MultiVector* dg_dx,
00155          Epetra_MultiVector* dg_dxdot,
00156          Epetra_MultiVector* dg_dxdotdot,
00157          Epetra_MultiVector* dg_dp);
00158 
00159 
00161     virtual void 
00162     postProcessResponses(const Epetra_Comm& comm, const Teuchos::RCP<Epetra_Vector>& g);
00163 
00165     virtual void 
00166     postProcessResponseDerivatives(const Epetra_Comm& comm, const Teuchos::RCP<Epetra_MultiVector>& gt);
00167 
00169     std::string getMode();
00170     bool pointIsInImagePtRegion(const double* p, double refZ) const;
00171     bool pointIsInAccumRegion(const double* p, double refZ) const;
00172     bool pointIsInLevelSetRegion(const double* p, double refZ) const;
00173     void addBeginPointData(const std::string& elementBlock, const double* p, double value);
00174     void addEndPointData(const std::string& elementBlock, const double* p, double value);
00175     void addImagePointData(const double* p, double value, double* grad);
00176     void addFinalImagePointData(const double* p, double value);
00177     void accumulatePointData(const double* p, double value, double* grad);
00178     void accumulateLevelSetData(const double* p, double value, double cellArea);
00179     double getSaddlePointWeight(const double* p) const;
00180     double getTotalSaddlePointWeight() const;
00181     const double* getSaddlePointPosition() const;
00182     double getCurrent(const double& lattTemp, const Teuchos::RCP<QCAD::MaterialDatabase>& materialDB) const;
00183     
00184   private:
00185 
00187     void initializeImagePoints(const double current_time, const Epetra_Vector* xdot,
00188              const Epetra_Vector& x, const Teuchos::Array<ParamVec>& p,
00189              Epetra_Vector& g, int dbMode);
00190     void initializeFinalImagePoints(const double current_time, const Epetra_Vector* xdot,
00191              const Epetra_Vector& x, const Teuchos::Array<ParamVec>& p,
00192              Epetra_Vector& g, int dbMode);
00193     void doNudgedElasticBand(const double current_time, const Epetra_Vector* xdot,
00194            const Epetra_Vector& x, const Teuchos::Array<ParamVec>& p,
00195            Epetra_Vector& g, int dbMode);
00196     void fillSaddlePointData(const double current_time, const Epetra_Vector* xdot,
00197            const Epetra_Vector& x, const Teuchos::Array<ParamVec>& p,
00198            Epetra_Vector& g, int dbMode);
00199 
00201     void doLevelSet(const double current_time,  const Epetra_Vector* xdot,
00202         const Epetra_Vector& x,  const Teuchos::Array<ParamVec>& p,
00203         Epetra_Vector& g, int dbMode);
00204     int FindSaddlePoint_LevelSet(std::vector<double>& allFieldVals,
00205            std::vector<double>* allCoords, std::vector<int>& ordering,
00206            double cutoffDistance, double cutoffFieldVal, double minDepth, int dbMode,
00207            Epetra_Vector& g);
00208 
00210     void getImagePointValues(const double current_time, const Epetra_Vector* xdot,
00211            const Epetra_Vector& x, const Teuchos::Array<ParamVec>& p,
00212            Epetra_Vector& g, double* globalPtValues, double* globalPtWeights,
00213            double* globalPtGrads, std::vector<mathVector> lastPositions, int dbMode);
00214     void getFinalImagePointValues(const double current_time, const Epetra_Vector* xdot,
00215            const Epetra_Vector& x, const Teuchos::Array<ParamVec>& p,
00216            Epetra_Vector& g, int dbMode);
00217     void writeOutput(int nIters);
00218     void initialIterationSetup(double& gradScale, double& springScale, int dbMode);
00219     void computeTangent(std::size_t i, mathVector& tangent, int dbMode);
00220     void computeClimbingForce(std::size_t i, const QCAD::mathVector& tangent, 
00221             const double& gradScale, QCAD::mathVector& force, int dbMode);
00222     void computeForce(std::size_t i, const QCAD::mathVector& tangent, 
00223           const std::vector<double>& springConstants,
00224           const double& gradScale,  const double& springScale, 
00225           QCAD::mathVector& force, double& dt, double& dt2, int dbMode);
00226 
00227     bool matchesCurrentResults(Epetra_Vector& g) const;
00228 
00229 
00231     SaddleValueResponseFunction(const SaddleValueResponseFunction&);
00232     
00234     SaddleValueResponseFunction& operator=(const SaddleValueResponseFunction&);
00235 
00237     double pointFn(double d, double radius) const;
00238 
00240     int getHighestPtIndex() const;
00241 
00243     Teuchos::RCP<const Epetra_Comm> comm;
00244 
00246     std::size_t numDims;
00247     std::size_t nImagePts;
00248     std::vector<nebImagePt> imagePts;
00249     std::vector<nebImagePt> finalPts;
00250     double imagePtSize;
00251     bool bClimbing, bAdaptivePointSize;
00252     double minAdaptivePointWt, maxAdaptivePointWt;
00253     double antiKinkFactor;
00254 
00256     bool bAggregateWorksets;
00257     std::vector<double> vFieldValues;
00258     std::vector<maxDimPt> vCoords;
00259     std::vector<maxDimPt> vGrads;
00260 
00262     std::vector<double> vlsFieldValues;
00263     std::vector<double> vlsCellAreas;
00264     std::vector<double> vlsCoords[MAX_DIMENSIONS];
00265     double fieldCutoffFctr;
00266     double minPoolDepthFctr;
00267     double distanceCutoffFctr;
00268     double levelSetRadius;
00269 
00270     // index into imagePts of the "found" saddle point
00271     int iSaddlePt;
00272     double returnFieldVal;  // value of the return field at the "found" saddle point
00273 
00274     double maxTimeStep, minTimeStep;
00275     double minSpringConstant, maxSpringConstant;
00276     std::size_t maxIterations;
00277     std::size_t backtraceAfterIters;
00278     double convergeTolerance;
00279 
00281     std::string beginRegionType, endRegionType; // "Point", "Element Block", or "Polygon"
00282     std::string beginElementBlock, endElementBlock;
00283     std::vector<mathVector> beginPolygon, endPolygon;
00284     bool saddleGuessGiven;
00285     mathVector saddlePointGuess;
00286     double shortenBeginPc, shortenEndPc;
00287 
00288     double zmin, zmax;  //defines lateral-volume region when numDims == 3
00289     double xmin, xmax, ymin, ymax; // dynamically adjusted box marking region containing image points
00290     bool bLockToPlane;
00291     double lockedZ;
00292 
00294     int maxFinalPts;
00295     
00296     double gfGridSpacing;  // grid spacing for GB-CBR calculation
00297     double fieldScaling;   // unscale the field specified by Field Name by Field Scaling Factor
00298     
00299     double initVds;        // initial Vds value in [V]
00300     double finalVds;       // final Vds value in [V]
00301     int stepsVds;          // number of steps from initial to final Vds
00302     bool bSweepVds;        // true if sweeping Vds, false if only using finalVds
00303     
00304     std::string gfEigensolver;  // Anasazi or tql2 eigensolver used for the GF-CBR calculation
00305 
00306     bool bGetCurrent;
00307     double current_Ecutoff_offset_from_Emax;
00308 
00310     mathVector imagePtValues;
00311     mathVector imagePtWeights;
00312     mathVector imagePtGradComps;
00313 
00314     mathVector finalPtValues;
00315     mathVector finalPtWeights;
00316 
00317 
00319     std::string mode;
00320 
00321     int  debugMode;
00322 
00323     std::string outputFilename;
00324     std::string debugFilename;
00325     int nEvery;
00326     bool appendOutput; // true => output is just appended if files exist
00327   };
00328 
00329   
00330 
00331 
00332 }
00333 
00334 #endif // QCAD_SADDLEVALUERESPONSEFUNCTION_HPP