Canoe
Comprehensive Atmosphere N' Ocean Engine
implicit_solver.hpp
Go to the documentation of this file.
1 #ifndef SRC_SNAP_IMPLICIT_IMPLICIT_SOLVER_HPP_
2 #define SRC_SNAP_IMPLICIT_IMPLICIT_SOLVER_HPP_
3 
4 // C/C++
5 #include <memory>
6 #include <string>
7 #include <vector>
8 
9 // Eigen
10 #include <Eigen/Core>
11 
12 // athena
13 #include <athena/athena.hpp>
14 // #include <bvals/bvals_interfaces.hpp>
15 
16 class MeshBlock;
17 class ParameterInput;
18 class BlockIndex;
19 class Coordinates;
20 class Mesh;
21 class BoundaryValues;
22 class EquationOfState;
23 class Thermodynamics;
24 
25 template <typename T>
26 class AthenaArray;
27 
29  friend class Hydro;
30 
31  public:
32  EIGEN_MAKE_ALIGNED_OPERATOR_NEW
33  // data
38  NeighborBlock tblock, bblock;
40 
41  // functions
42  ImplicitSolver(MeshBlock *pmb, ParameterInput *pin);
44  void SetDirection(CoordinateDirection dir);
45  int GetImplicitFlag() { return implicit_flag; }
46 
47  // utility functions
48  void FindNeighbors();
49  int CreateMPITag(int lid, int bufid, std::string phy);
50 
51  // correction methods
53  Real dt);
55  Real dt);
56 
57  // tri-diagonal solver
58  template <typename T1, typename T2>
59  void ForwardSweep(std::vector<T1> &a, std::vector<T1> &b, std::vector<T1> &c,
60  std::vector<T2> &delta, std::vector<T2> &corr, Real dt,
61  int k, int j, int il, int iu);
62 
63  template <typename T1, typename T2>
64  void BackwardSubstitution(std::vector<T1> &a, std::vector<T2> &delta, int kl,
65  int ku, int jl, int ju, int il, int iu);
66 
67  // periodic solver
68  template <typename T1, typename T2>
69  void PeriodicForwardSweep(std::vector<T1> &a, std::vector<T1> &b,
70  std::vector<T1> &c, std::vector<T2> &corr, Real dt,
71  int k, int j, int il, int iu);
72 
73  template <typename T1, typename T2>
74  void PeriodicBackwardSubstitution(std::vector<T1> &a, std::vector<T1> &c,
75  std::vector<T2> &delta, int kl, int ku,
76  int jl, int ju, int il, int iu);
77 
78  // forcing jacobians
79  template <typename T>
80  void JacobianGravityCoriolis(T &jac, Real const prim[], int k, int j, int i);
81 
82  // communications
83  /*void SynchronizeConserved(AthenaArray<Real> const& du,
84  int kl, int ku, int jl, int ju, int is, int ie);
85  void WaitToFinishSync(int kl, int ku, int jl, int ju, int is, int ie);*/
86 
87  template <typename T>
88  void SendBuffer(T const &a, int k, int j, NeighborBlock nb);
89 
90  template <typename T1, typename T2>
91  void SendBuffer(T1 const &a, T2 const &b, int k, int j, NeighborBlock nb);
92 
93  template <typename T1, typename T2, typename T3, typename T4, typename T5,
94  typename T6>
95  void SendBuffer(T1 const &a, T2 const &b, T3 const &c, T4 const &d,
96  T5 const &e, T6 const &f, int k, int j, NeighborBlock ntop);
97 
98  template <typename T1, typename T2, typename T3, typename T4, typename T5,
99  typename T6, typename T7>
100  void SendBuffer(T1 const &a, T2 const &b, T3 const &c, T4 const &d,
101  T5 const &e, T6 const &f, T7 const &g, int k, int j,
102  NeighborBlock ntop);
103 
104  template <typename T>
105  void RecvBuffer(T &a, int k, int j, NeighborBlock nb);
106 
107  template <typename T1, typename T2>
108  void RecvBuffer(T1 &a, T2 &b, int k, int j, NeighborBlock nb);
109 
110  template <typename T1, typename T2, typename T3, typename T4, typename T5,
111  typename T6>
112  void RecvBuffer(T1 &a, T2 &b, T3 &c, T4 &d, T5 &e, T6 &f, int k, int j,
113  NeighborBlock nb);
114 
115  template <typename T1, typename T2, typename T3, typename T4, typename T5,
116  typename T6, typename T7>
117  void RecvBuffer(T1 &a, T2 &b, T3 &c, T4 &d, T5 &e, T6 &f, T7 &g, int k, int j,
118  NeighborBlock nb);
119 
120  template <typename T1, typename T2>
121  void SaveCoefficients(std::vector<T1> &a, std::vector<T2> &b, int k, int j,
122  int il, int iu);
123 
124  template <typename T1, typename T2, typename T3, typename T4>
125  void SaveCoefficients(std::vector<T1> &a, std::vector<T2> &b,
126  std::vector<T3> &c, std::vector<T4> &d, int k, int j,
127  int il, int iu);
128 
129  template <typename T1, typename T2>
130  void LoadCoefficients(std::vector<T1> &a, std::vector<T2> &b, int k, int j,
131  int il, int iu);
132 
133  template <typename T1, typename T2, typename T3, typename T4>
134  void LoadCoefficients(std::vector<T1> &a, std::vector<T2> &b,
135  std::vector<T3> &c, std::vector<T4> &d, int k, int j,
136  int il, int iu);
137 
138  // template<typename T>
139  // void LoadForcingJacobian(T &phi, int k, int j ,int i, CoordinateDirection
140  // dir);
141 
142  private:
143  MeshBlock const *pmy_block_;
144 
145  CoordinateDirection mydir_;
146  // Real *usend_top_, *urecv_top_; // MPI data buffer
147  // Real *usend_bot_, *urecv_bot_; // MPI data buffer
148 
149  Real ***buffer_; // MPI data buffer
150  // Real ****jacobian_; // archive of forcing jacobian
151  AthenaArray<Real> du_; // stores implicit solution
153  coefficients_; // archive of coefficients in the tri-diagonal matrix
154 
155  Eigen::Matrix<Real, 5, 5> p2_, p3_; // perturbation matrices
156 
157 #ifdef MPI_PARALLEL
158  MPI_Request **req_send_data1_;
159  MPI_Request **req_send_data2_;
160  MPI_Request **req_send_data6_;
161  MPI_Request **req_send_data7_;
162  // MPI_Request req_send_sync_top_;
163  // MPI_Request req_send_sync_bot_;
164  // MPI_Request req_recv_sync_top_;
165  // MPI_Request req_recv_sync_bot_;
166 #endif
167 };
168 
169 using ImplicitSolverPtr = std::shared_ptr<ImplicitSolver>;
170 
171 #endif // SRC_SNAP_IMPLICIT_IMPLICIT_SOLVER_HPP_
EIGEN_MAKE_ALIGNED_OPERATOR_NEW bool has_top_neighbor
void FullCorrection(AthenaArray< Real > &du, AthenaArray< Real > const &w, Real dt)
int CreateMPITag(int lid, int bufid, std::string phy)
void ForwardSweep(std::vector< T1 > &a, std::vector< T1 > &b, std::vector< T1 > &c, std::vector< T2 > &delta, std::vector< T2 > &corr, Real dt, int k, int j, int il, int iu)
friend class Hydro
NeighborBlock bblock
EIGEN_MAKE_ALIGNED_OPERATOR_NEW bool has_bot_neighbor
void RecvBuffer(T &a, int k, int j, NeighborBlock nb)
NeighborBlock tblock
Eigen::Matrix< Real, 5, 5 > p3_
void PartialCorrection(AthenaArray< Real > &du, AthenaArray< Real > const &w, Real dt)
CoordinateDirection mydir_
AthenaArray< Real > coefficients_
void JacobianGravityCoriolis(T &jac, Real const prim[], int k, int j, int i)
void PeriodicBackwardSubstitution(std::vector< T1 > &a, std::vector< T1 > &c, std::vector< T2 > &delta, int kl, int ku, int jl, int ju, int il, int iu)
void BackwardSubstitution(std::vector< T1 > &a, std::vector< T2 > &delta, int kl, int ku, int jl, int ju, int il, int iu)
void SetDirection(CoordinateDirection dir)
void LoadCoefficients(std::vector< T1 > &a, std::vector< T2 > &b, int k, int j, int il, int iu)
AthenaArray< Real > du_
void PeriodicForwardSweep(std::vector< T1 > &a, std::vector< T1 > &b, std::vector< T1 > &c, std::vector< T2 > &corr, Real dt, int k, int j, int il, int iu)
void SaveCoefficients(std::vector< T1 > &a, std::vector< T2 > &b, int k, int j, int il, int iu)
Eigen::Matrix< Real, 5, 5 > p2_
MeshBlock const * pmy_block_
ImplicitSolver(MeshBlock *pmb, ParameterInput *pin)
void SendBuffer(T const &a, int k, int j, NeighborBlock nb)
std::shared_ptr< ImplicitSolver > ImplicitSolverPtr