1 // Copyright (c) 2017-2021, Lawrence Livermore National Security, LLC and
2 // other Axom Project Developers. See the top-level LICENSE file for details.
3 //
4 // SPDX-License-Identifier: (BSD-3-Clause)
5 
6 
7 // Hydro driver
8 // Fri Mar 27 14:14:28 PDT 2015
9 
10 // HYDRO C: Physics-wise, this hydro is still pretty dumb: multiple
11 // materials, gamma-law gas.  The mesh is an arbitrary collection of
12 // polygons.  Architecturally, we use as few temporaries as possible.
13 // No sweet returning of objects when you can instead pass in an array
14 // and modify it in place.  No STL containers.  Explicit memory
15 // allocation of all the memory we ever need at problem start, and
16 // explicit deletes when the hydro object is destroyed.
17 
18 #include "State.hpp"
19 
20 #include "TinyHydroTypes.hpp"
21 
22 namespace tinyHydro {
23 
24   class Hydro
25   {
26   public:
27 
28   public:
29     Hydro(State * s); // using a pointer here with PyBindGen prevents
30                       // State being garbage collected by Python
31     ~Hydro();
32 
33     // problem control and advance methods
34     void      step(double dt);
35     double    newDT();
36     void      steps(int numSteps);
37     void      advance(double stopTime);
38     void      initialize(void); // call once before taking steps
39     void      setBC(const char * boundary, double xVel, double yVel);
40 
41     // compute time derivatives from first arg, store rates in 2nd arg
42     void      calcDerivs(const State & state, State & dState, double dt);
43 
44     // physics methods
45     void      calcQ(const State & state);
46     void      calcForce(const State & state);
47     void      calcDivU(const NodalVectorField& u);
48     void      applyVelocityBC(NodalVectorField& u);
49     double    totalEnergy(const State & s) const;
50     double    totalEnergy(void) const;
51     void      calcMaxCs(void);
52 
53     // update all the algebraically determinable stuff, like density
54     void      updateConstitutives(State & state);
55 
56     // setter methods (used from Python)
setState(const State s)57     void      setState(const State s) {state = s; }
58 
59     // getter methods
getState(void)60     State *   getState(void) {return &state; }
getQ(int i)61     double    getQ(int i) {return Q[i]; }
getForce(int i)62     VectorXY  getForce(int i) {return force[i]; }
63     int       numBCnodes(int bc);
64     int       bcNode(int bc, int node);
65 
fuzzyEqual(double a,double b,double tol=1.0e-15)66     bool      fuzzyEqual(double a, double b, double tol = 1.0e-15)
67     {
68       return (fabs(a - b) < tol);
69     }
70 
71   public:
72     PolygonMeshXY&       mesh;
73 
74     State&               state;
75     State halfStep;                         // re-used temporary for RK2
76     State dState;                           // re-used temporary for RK2
77 
78     ZonalScalarField Q;                     // re-used for Q
79     ZonalScalarField divu;                  // for dt, Q, and some PdV work schemes
80 
81     NodalVectorField force;                 // re-used for force calc
82     NodalVectorField halfStepVelocity;      // re-used for work calc
83 
84     BoundaryEdgeSet boundaryEdgeSet;                                  // SLAM -- new set for domain boundaries (of size 4)
85     BoundaryEdgeVectorField bcVelocity;     // velocity BC values     // SLAM -- map to velocity boundary conditions for each edge
86     NodeSubset bcNodes[NUM_DOMAIN_BOUNDARIES];                        // SLAM -- indirection sets of nodes on each boundary edge
87 
88 
89     ZonalScalarField*    initialMass;   // per material initial mass of zones
90 
91     ZonalScalarField totalMass;         // overall mass in zone
92     NodalScalarField nodeMass;          // also keep node masses
93 
94     ZonalScalarField*    pressure;      // per material zone pressure
95 
96     ZonalScalarField totalPressure;     // sum of pressures in a zone
97     ZonalScalarField maxCs;             //
98 
99 
100     int cycle;
101     int dtZone;
102     double time;
103     double Cq, Cl;                      // Q coefficients
104     double cfl;
105     int printCycle;
106 
107   private:
108     // Unused variables...
109     // int numBCs;                                                      // <UNUSED>?
110 
111 #ifdef TINY_HYDRO_USE_CORNER_DATA
112     VectorXY * cornerforce; // re-used for force calc                // <UNUSED>?
113 #endif
114   };
115 
116 } // end namespace tinyHydro
117