diff --git a/algebra/TinyMatrix.hpp b/algebra/TinyMatrix.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..94bf1ad7d020913999f01e88af84e39e6dd5f97c
--- /dev/null
+++ b/algebra/TinyMatrix.hpp
@@ -0,0 +1,249 @@
+#ifndef TINY_MATRIX_HPP
+#define TINY_MATRIX_HPP
+
+#include <cassert>
+#include <iostream>
+
+#include <Types.hpp>
+#include <TinyVector.hpp>
+
+template <size_t N, typename T=double>
+class TinyMatrix
+{
+private:
+ T m_values[N*N];
+ static_assert((N>0),"TinyMatrix size must be strictly positive");
+
+ KOKKOS_INLINE_FUNCTION
+ size_t _index(const size_t& i, const size_t& j) const
+ {
+ return std::move(i*N+j);
+ }
+
+ void _unpackVariadicInput(const T& t)
+ {
+ m_values[N*N-1] = t;
+ }
+
+ template <typename... Args>
+ void _unpackVariadicInput(const T& t, Args&&... args)
+ {
+ m_values[N*N-1-sizeof...(args)] = t;
+ this->_unpackVariadicInput(args...);
+ }
+
+public:
+
+ KOKKOS_INLINE_FUNCTION
+ friend TinyMatrix operator*(const T& t, const TinyMatrix& A)
+ {
+ TinyMatrix tA;
+ for (size_t i=0; i<N*N; ++i) {
+ tA.m_values[i] = t * A.m_values[i];
+ }
+ return std::move(tA);
+ }
+
+
+ KOKKOS_INLINE_FUNCTION
+ TinyMatrix operator*(const TinyMatrix& B) const
+ {
+ const TinyMatrix& A = *this;
+ TinyMatrix AB;
+ for (size_t i=0; i<N; ++i) {
+ for (size_t j=0; j<N; ++j) {
+ T sum = A(i,0)*B(0,j);
+ for (size_t k=1; k<N; ++k) {
+ sum += A(i,k)*B(k,j);
+ }
+ AB(i,j) = sum;
+ }
+ }
+ return std::move(AB);
+ }
+
+
+ KOKKOS_INLINE_FUNCTION
+ TinyVector<N,T> operator*(const TinyVector<N,T>& x) const
+ {
+ const TinyMatrix& A = *this;
+ TinyVector<N,T> Ax;
+ for (size_t i=0; i<N; ++i) {
+ T sum = A(i,0)*x[0];
+ for (size_t j=1; j<N; ++j) {
+ sum += A(i,j)*x[j];
+ }
+ Ax[i] = sum;
+ }
+ return std::move(Ax);
+ }
+
+
+ KOKKOS_INLINE_FUNCTION
+ friend std::ostream& operator<<(std::ostream& os, const TinyMatrix& A)
+ {
+ os << '[';
+ for (size_t i=0; i<N; ++i) {
+ os << '(' << A(i,0);
+ for (size_t j=1; j<N; ++j) {
+ os << ',' << A(i,j);
+ }
+ os << ')';
+ }
+ os << ']';
+ return os;
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ TinyMatrix operator+(const TinyMatrix& A) const
+ {
+ TinyMatrix sum;
+ for (size_t i=0; i<N*N; ++i) {
+ sum.m_values[i] = m_values[i]+A.m_values[i];
+ }
+ return std::move(sum);
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ TinyMatrix operator-(const TinyMatrix& A) const
+ {
+ TinyMatrix difference;
+ for (size_t i=0; i<N*N; ++i) {
+ difference.m_values[i] = m_values[i]-A.m_values[i];
+ }
+ return std::move(difference);
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ TinyMatrix& operator+=(const TinyMatrix& A)
+ {
+ for (size_t i=0; i<N*N; ++i) {
+ m_values[i] += A.m_values[i];
+ }
+ return *this;
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ TinyMatrix& operator-=(const TinyMatrix& A)
+ {
+ for (size_t i=0; i<N*N; ++i) {
+ m_values[i] -= A.m_values[i];
+ }
+ return *this;
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ T& operator()(const size_t& i, const size_t& j)
+ {
+ assert((i<N) and (j<N));
+ return m_values[_index(i,j)];
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ const T& operator()(const size_t& i, const size_t& j) const
+ {
+ assert((i<N) and (j<N));
+ return m_values[_index(i,j)];
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ TinyMatrix& operator=(const ZeroType& z)
+ {
+ static_assert(std::is_arithmetic<T>(),"Cannot assign 'zero' value for non-arithmetic types");
+ for (size_t i=0; i<N*N; ++i) {
+ m_values[i] = 0;
+ }
+ return *this;
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ TinyMatrix& operator=(const IdentityType& I)
+ {
+ static_assert(std::is_arithmetic<T>(),"Cannot assign 'identity' value for non-arithmetic types");
+ for (size_t i=0; i<N; ++i) {
+ for (size_t j=0; j<N; ++j) {
+ m_values[_index(i,j)] = (i==j) ? 1 : 0;
+ }
+ }
+ return *this;
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ const TinyMatrix& operator=(const TinyMatrix& A)
+ {
+ for (size_t i=0; i<N*N; ++i) {
+ m_values[i] = A.m_values[i];
+ }
+ return *this;
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ TinyMatrix& operator=(TinyMatrix&& A) = default;
+
+ template <typename... Args>
+ KOKKOS_INLINE_FUNCTION
+ TinyMatrix(const T& t, Args&&... args)
+ {
+ static_assert(sizeof...(args)==N*N-1, "wrong number of parameters");
+ this->_unpackVariadicInput(t, args...);
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ TinyMatrix()
+ {
+ ;
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ TinyMatrix(const ZeroType& z)
+ {
+ static_assert(std::is_arithmetic<T>(),"Cannot construct from 'zero' value for non-arithmetic types");
+ for (size_t i=0; i<N*N; ++i) {
+ m_values[i] = 0;
+ }
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ TinyMatrix(const IdentityType& I)
+ {
+ static_assert(std::is_arithmetic<T>(),"Cannot construct from 'identity' value for non-arithmetic types");
+ for (size_t i=0; i<N; ++i) {
+ for (size_t j=0; j<N; ++j) {
+ m_values[_index(i,j)] = (i==j) ? 1 : 0;
+ }
+ }
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ TinyMatrix(const TinyMatrix& A)
+ {
+ for (size_t i=0; i<N*N; ++i) {
+ m_values[i] = A.m_values[i];
+ }
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ TinyMatrix(TinyMatrix&& A) = default;
+
+ KOKKOS_INLINE_FUNCTION
+ ~TinyMatrix()
+ {
+ ;
+ }
+};
+
+template <size_t N, typename T>
+KOKKOS_INLINE_FUNCTION
+TinyMatrix<N,T> tensorProduct(const TinyVector<N,T>& x,
+ const TinyVector<N,T>& y)
+{
+ TinyMatrix<N,T> A;
+ for (size_t i=0; i<N; ++i) {
+ for (size_t j=0; j<N; ++j) {
+ A(i,j) = x[i]*y[j];
+ }
+ }
+ return std::move(A);
+}
+
+#endif // TINYMATRIX_HPP
diff --git a/algebra/TinyVector.hpp b/algebra/TinyVector.hpp
index 95dfda71257c53d514ce2f5dd89c647b6c382f4b..9b2ca1619f0de8e2b811d41ae7be7f711c68a436 100644
--- a/algebra/TinyVector.hpp
+++ b/algebra/TinyVector.hpp
@@ -27,7 +27,7 @@ private:
public:
KOKKOS_INLINE_FUNCTION
- T operator,(const TinyVector& v)
+ T operator,(const TinyVector& v) const
{
T t = m_values[0]*v.m_values[0];
for (size_t i=1; i<N; ++i) {
diff --git a/experimental/AcousticSolverTest.cpp b/experimental/AcousticSolverTest.cpp
index fd5e9f994ffcc6c9cbfb65614632beed246c2d4f..67a54b78881ab3d65f1c5767019eac4dced40fce 100644
--- a/experimental/AcousticSolverTest.cpp
+++ b/experimental/AcousticSolverTest.cpp
@@ -53,13 +53,13 @@ AcousticSolverTest::computeRhoCj(const Kokkos::View<const double*>& rhoj,
}
KOKKOS_INLINE_FUNCTION
-const Kokkos::View<const double*[2]>
+const Kokkos::View<const Rdd*[2]>
AcousticSolverTest::computeAjr(const Kokkos::View<const double*>& rhocj,
- const Kokkos::View<const double*[2]>& Cjr)
+ const Kokkos::View<const Rd*[2]>& Cjr)
{
Kokkos::parallel_for(m_nj, KOKKOS_LAMBDA(const int& j) {
for (int r=0; r<2; ++r) {
- m_Ajr(j,r) = rhocj(j)*Cjr(j,r)*Cjr(j,r);
+ m_Ajr(j,r) = tensorProduct(rhocj(j)*Cjr(j,r), Cjr(j,r));
}
});
@@ -67,14 +67,14 @@ AcousticSolverTest::computeAjr(const Kokkos::View<const double*>& rhocj,
}
KOKKOS_INLINE_FUNCTION
-const Kokkos::View<const double*>
-AcousticSolverTest::computeAr(const Kokkos::View<const double*[2]>& Ajr,
+const Kokkos::View<const Rdd*>
+AcousticSolverTest::computeAr(const Kokkos::View<const Rdd*[2]>& Ajr,
const Kokkos::View<const int*[2]>& node_cells,
const Kokkos::View<const int*[2]>& node_cell_local_node,
const Kokkos::View<const int*>& node_nb_cells)
{
Kokkos::parallel_for(m_nr, KOKKOS_LAMBDA(const int& r) {
- double sum = 0;
+ Rdd sum = zero;
for (int j=0; j<node_nb_cells(r); ++j) {
const int J = node_cells(r,j);
const int R = node_cell_local_node(r,j);
@@ -87,56 +87,56 @@ AcousticSolverTest::computeAr(const Kokkos::View<const double*[2]>& Ajr,
}
KOKKOS_INLINE_FUNCTION
-const Kokkos::View<const double*>
-AcousticSolverTest::computeBr(const Kokkos::View<const double*[2]>& Ajr,
- const Kokkos::View<const double*[2]>& Cjr,
- const Kokkos::View<const double*>& uj,
+const Kokkos::View<const Rd*>
+AcousticSolverTest::computeBr(const Kokkos::View<const Rdd*[2]>& Ajr,
+ const Kokkos::View<const Rd*[2]>& Cjr,
+ const Kokkos::View<const Rd*>& uj,
const Kokkos::View<const double*>& pj,
const Kokkos::View<const int*[2]>& node_cells,
const Kokkos::View<const int*[2]>& node_cell_local_node,
const Kokkos::View<const int*>& node_nb_cells)
{
Kokkos::parallel_for(m_nr, KOKKOS_LAMBDA(const int& r) {
- double sum = 0;
+ Rd& br = m_br(r);
+ br = zero;
for (int j=0; j<node_nb_cells(r); ++j) {
const int J = node_cells(r,j);
const int R = node_cell_local_node(r,j);
- sum += Ajr(J,R)*uj(J) + Cjr(J,R)*pj(J);
+ br += Ajr(J,R)*uj(J) + pj(J)*Cjr(J,R);
}
- m_br(r) = sum;
});
return m_br;
}
KOKKOS_INLINE_FUNCTION
-Kokkos::View<double*>
-AcousticSolverTest::computeUr(const Kokkos::View<const double*>& Ar,
- const Kokkos::View<const double*>& br)
+Kokkos::View<Rd*>
+AcousticSolverTest::computeUr(const Kokkos::View<const Rdd*>& Ar,
+ const Kokkos::View<const Rd*>& br)
{
inverse(Ar, m_inv_Ar);
- const Kokkos::View<const double*> invAr = m_inv_Ar;
+ const Kokkos::View<const Rdd*> invAr = m_inv_Ar;
Kokkos::parallel_for(m_nr, KOKKOS_LAMBDA(const int& r) {
- m_ur[r]=br(r)*invAr(r);
+ m_ur[r]=invAr(r)*br(r);
});
- m_ur[0]=0;
- m_ur[m_nr-1]=0;
+ m_ur[0]=zero;
+ m_ur[m_nr-1]=zero;
return m_ur;
}
KOKKOS_INLINE_FUNCTION
-Kokkos::View<double*[2]>
-AcousticSolverTest::computeFjr(const Kokkos::View<const double*[2]>& Ajr,
- const Kokkos::View<const double*>& ur,
- const Kokkos::View<const double*[2]>& Cjr,
- const Kokkos::View<const double*>& uj,
+Kokkos::View<Rd*[2]>
+AcousticSolverTest::computeFjr(const Kokkos::View<const Rdd*[2]>& Ajr,
+ const Kokkos::View<const Rd*>& ur,
+ const Kokkos::View<const Rd*[2]>& Cjr,
+ const Kokkos::View<const Rd*>& uj,
const Kokkos::View<const double*>& pj,
const Kokkos::View<const int*[2]>& cell_nodes)
{
Kokkos::parallel_for(m_nj, KOKKOS_LAMBDA(const int& j) {
for (int r=0; r<2; ++r) {
- m_Fjr(j,r) = Ajr(j,r)*(uj(j)-ur(cell_nodes(j,r)))+Cjr(j,r)*pj(j);
+ m_Fjr(j,r) = Ajr(j,r)*(uj(j)-ur(cell_nodes(j,r)))+pj(j)*Cjr(j,r);
}
});
@@ -169,26 +169,37 @@ AcousticSolverTest::inverse(const Kokkos::View<const double*>& x,
}
KOKKOS_INLINE_FUNCTION
-void AcousticSolverTest::computeExplicitFluxes(const Kokkos::View<const double*>& xr,
- const Kokkos::View<const double*>& xj,
+void
+AcousticSolverTest::inverse(const Kokkos::View<const Rdd*>& A,
+ Kokkos::View<Rdd*>& inv_A) const
+{
+ Kokkos::parallel_for(A.size(), KOKKOS_LAMBDA(const int& r) {
+ inv_A(r) = Rdd{1./(A(r)(0,0))};
+ });
+}
+
+
+KOKKOS_INLINE_FUNCTION
+void AcousticSolverTest::computeExplicitFluxes(const Kokkos::View<const Rd*>& xr,
+ const Kokkos::View<const Rd*>& xj,
const Kokkos::View<const double*>& rhoj,
- const Kokkos::View<const double*>& uj,
+ const Kokkos::View<const Rd*>& uj,
const Kokkos::View<const double*>& pj,
const Kokkos::View<const double*>& cj,
const Kokkos::View<const double*>& Vj,
- const Kokkos::View<const double*[2]>& Cjr,
+ const Kokkos::View<const Rd*[2]>& Cjr,
const Kokkos::View<const int*[2]>& cell_nodes,
const Kokkos::View<const int*[2]>& node_cells,
const Kokkos::View<const int*>& node_nb_cells,
const Kokkos::View<const int*[2]>& node_cell_local_node,
- Kokkos::View<double*>& ur,
- Kokkos::View<double*[2]>& Fjr)
+ Kokkos::View<Rd*>& ur,
+ Kokkos::View<Rd*[2]>& Fjr)
{
const Kokkos::View<const double*> rhocj = computeRhoCj(rhoj, cj);
- const Kokkos::View<const double*[2]> Ajr = computeAjr(rhocj, Cjr);
+ const Kokkos::View<const Rdd*[2]> Ajr = computeAjr(rhocj, Cjr);
- const Kokkos::View<const double*> Ar = computeAr(Ajr, node_cells, node_cell_local_node, node_nb_cells);
- const Kokkos::View<const double*> br = computeBr(Ajr, Cjr, uj, pj,
+ const Kokkos::View<const Rdd*> Ar = computeAr(Ajr, node_cells, node_cell_local_node, node_nb_cells);
+ const Kokkos::View<const Rd*> br = computeBr(Ajr, Cjr, uj, pj,
node_cells, node_cell_local_node, node_nb_cells);
ur = computeUr(Ar, br);
@@ -207,10 +218,10 @@ AcousticSolverTest::AcousticSolverTest(const long int& nj)
m_rhocj("rho_c", m_nj),
m_Vj_over_cj("Vj_over_cj", m_nj)
{
- Kokkos::View<double*> xj("xj",m_nj);
+ Kokkos::View<Rd*> xj("xj",m_nj);
Kokkos::View<double*> rhoj("rhoj",m_nj);
- Kokkos::View<double*> uj("uj",m_nj);
+ Kokkos::View<Rd*> uj("uj",m_nj);
Kokkos::View<double*> Ej("Ej",m_nj);
Kokkos::View<double*> ej("ej",m_nj);
@@ -220,12 +231,12 @@ AcousticSolverTest::AcousticSolverTest(const long int& nj)
Kokkos::View<double*> cj("cj",m_nj);
Kokkos::View<double*> mj("mj",m_nj);
- Kokkos::View<double*> xr("xr", m_nr);
+ Kokkos::View<Rd*> xr("xr", m_nr);
const double delta_x = 1./m_nj;
Kokkos::parallel_for(m_nr, KOKKOS_LAMBDA(const int& r){
- xr[r] = r*delta_x;
+ xr[r][0] = r*delta_x;
});
Kokkos::View<int*[2]> cell_nodes("cell_nodes",m_nj,2);
@@ -266,12 +277,19 @@ AcousticSolverTest::AcousticSolverTest(const long int& nj)
xj[j] = 0.5*(xr[cell_nodes(j,0)]+xr[cell_nodes(j,1)]);
});
+ Kokkos::View<Rd*[2]> Cjr("Cjr",m_nj);
+ Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j) {
+ Cjr(j,0)=-1;
+ Cjr(j,1)= 1;
+ });
+
Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- Vj[j] = xr[cell_nodes(j,1)]-xr[cell_nodes(j,0)];
+ Vj[j] = (xr[cell_nodes(j,1)], Cjr(j,1))
+ + (xr[cell_nodes(j,0)], Cjr(j,0));
});
Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- if (xj[j]<0.5) {
+ if (xj[j][0]<0.5) {
rhoj[j]=1;
} else {
rhoj[j]=0.125;
@@ -279,13 +297,17 @@ AcousticSolverTest::AcousticSolverTest(const long int& nj)
});
Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- if (xj[j]<0.5) {
+ if (xj[j][0]<0.5) {
pj[j]=1;
} else {
pj[j]=0.1;
}
});
+ Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
+ uj[j] = zero;
+ });
+
Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
gammaj[j] = 1.4;
});
@@ -295,7 +317,7 @@ AcousticSolverTest::AcousticSolverTest(const long int& nj)
block_eos.updateEandCFromRhoP();
Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- Ej[j] = ej[j]+0.5*uj[j]*uj[j];
+ Ej[j] = ej[j]+0.5*(uj[j],uj[j]);
});
Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
@@ -311,24 +333,6 @@ AcousticSolverTest::AcousticSolverTest(const long int& nj)
int itermax=std::numeric_limits<int>::max();
int iteration=0;
- Kokkos::View<TinyVector<2>*[2]> Bjr("Cjr",m_nj);
- Bjr(0,0)=zero;
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j) {
- Bjr(j,0)[0]=j;
- Bjr(j,0)[1]=2;
-
- Bjr(j,1)=zero;
-});
-
- for (size_t j=0; j<m_nj; ++j) {
- std::cout << Bjr(j,0) << " | " << Bjr(j,1) << '\n';
- }
- Kokkos::View<double*[2]> Cjr("Cjr",m_nj);
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j) {
- Cjr(j,0)=-1;
- Cjr(j,1)= 1;
- });
-
while((t<tmax) and (iteration<itermax)) {
double dt = 0.4*acoustic_dt(Vj, cj);
if (t+dt<tmax) {
@@ -343,23 +347,23 @@ AcousticSolverTest::AcousticSolverTest(const long int& nj)
cell_nodes, node_cells, node_nb_cells, node_cell_local_node,
m_ur, m_Fjr);
- const Kokkos::View<const double*[2]> Fjr = m_Fjr;
- const Kokkos::View<const double*> ur = m_ur;
+ const Kokkos::View<const Rd*[2]> Fjr = m_Fjr;
+ const Kokkos::View<const Rd*> ur = m_ur;
Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j) {
- double momentum_fluxes = 0;
+ Rd momentum_fluxes = zero;
double energy_fluxes = 0;
for (int R=0; R<2; ++R) {
const int r=cell_nodes(j,R);
- momentum_fluxes += Fjr(j,R);
- energy_fluxes += Fjr(j,R)*ur[r];
+ momentum_fluxes += Fjr(j,R);
+ energy_fluxes += ((Fjr(j,R), ur[r]));
}
- uj[j] -= dt*inv_mj[j]*(momentum_fluxes);
- Ej[j] -= dt*inv_mj[j]*(energy_fluxes);
+ uj[j] -= (dt*inv_mj[j]) * momentum_fluxes;
+ Ej[j] -= (dt*inv_mj[j]) * energy_fluxes;
});
Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j) {
- ej[j] = Ej[j] - 0.5 * uj[j]*uj[j];
+ ej[j] = Ej[j] - 0.5 * (uj[j],uj[j]);
});
Kokkos::parallel_for(m_nr, KOKKOS_LAMBDA(const int& r){
@@ -368,9 +372,13 @@ AcousticSolverTest::AcousticSolverTest(const long int& nj)
Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
xj[j] = 0.5*(xr[j]+xr[j+1]);
- Vj[j] = xr[j+1]-xr[j];
});
+ Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
+ Vj[j] = (xr[cell_nodes(j,1)], Cjr(j,1))
+ + (xr[cell_nodes(j,0)], Cjr(j,0));
+ });
+
Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
rhoj[j] = mj[j]/Vj[j];
});
@@ -383,7 +391,7 @@ AcousticSolverTest::AcousticSolverTest(const long int& nj)
// {
// std::ofstream fout("rho");
// for (int j=0; j<nj; ++j) {
- // fout << xj[j] << ' ' << rhoj[j] << '\n';
+ // fout << xj[j][0] << ' ' << rhoj[j] << '\n';
// }
// }
diff --git a/experimental/AcousticSolverTest.hpp b/experimental/AcousticSolverTest.hpp
index cd6ba4b64dd1372d07095a306aa2ccb40e03a637..209a18747bb42eb467e6a992faaf8f64e706f856 100644
--- a/experimental/AcousticSolverTest.hpp
+++ b/experimental/AcousticSolverTest.hpp
@@ -4,79 +4,87 @@
#include <Kokkos_Core.hpp>
#include <TinyVector.hpp>
+#include <TinyMatrix.hpp>
+
+constexpr static size_t dimension = 1;
+
+typedef TinyVector<dimension> Rd;
+typedef TinyMatrix<dimension> Rdd;
class AcousticSolverTest
{
private:
- constexpr static size_t dimension = 1;
-
+
inline const Kokkos::View<const double*>
computeRhoCj(const Kokkos::View<const double*>& rhoj,
const Kokkos::View<const double*>& cj);
- inline const Kokkos::View<const double*[2]>
+ inline const Kokkos::View<const Rdd*[2]>
computeAjr(const Kokkos::View<const double*>& rhocj,
- const Kokkos::View<const double*[2]>& Cjr);
+ const Kokkos::View<const Rd*[2]>& Cjr);
- inline const Kokkos::View<const double*>
- computeAr(const Kokkos::View<const double*[2]>& Ajr,
+ inline const Kokkos::View<const Rdd*>
+ computeAr(const Kokkos::View<const Rdd*[2]>& Ajr,
const Kokkos::View<const int*[2]>& node_cells,
const Kokkos::View<const int*[2]>& node_cell_local_node,
const Kokkos::View<const int*>& node_nb_cells);
- inline const Kokkos::View<const double*>
- computeBr(const Kokkos::View<const double*[2]>& Ajr,
- const Kokkos::View<const double*[2]>& Cjr,
- const Kokkos::View<const double*>& uj,
+ inline const Kokkos::View<const Rd*>
+ computeBr(const Kokkos::View<const Rdd*[2]>& Ajr,
+ const Kokkos::View<const Rd*[2]>& Cjr,
+ const Kokkos::View<const Rd*>& uj,
const Kokkos::View<const double*>& pj,
const Kokkos::View<const int*[2]>& node_cells,
const Kokkos::View<const int*[2]>& node_cell_local_node,
const Kokkos::View<const int*>& node_nb_cells);
- Kokkos::View<double*>
- computeUr(const Kokkos::View<const double*>& Ar,
- const Kokkos::View<const double*>& br);
+ Kokkos::View<Rd*>
+ computeUr(const Kokkos::View<const Rdd*>& Ar,
+ const Kokkos::View<const Rd*>& br);
- Kokkos::View<double*[2]>
- computeFjr(const Kokkos::View<const double*[2]>& Ajr,
- const Kokkos::View<const double*>& ur,
- const Kokkos::View<const double*[2]>& Cjr,
- const Kokkos::View<const double*>& uj,
+ Kokkos::View<Rd*[2]>
+ computeFjr(const Kokkos::View<const Rdd*[2]>& Ajr,
+ const Kokkos::View<const Rd*>& ur,
+ const Kokkos::View<const Rd*[2]>& Cjr,
+ const Kokkos::View<const Rd*>& uj,
const Kokkos::View<const double*>& pj,
const Kokkos::View<const int*[2]>& cell_nodes);
+ void inverse(const Kokkos::View<const Rdd*>& A,
+ Kokkos::View<Rdd*>& inv_A) const;
+
void inverse(const Kokkos::View<const double*>& x,
Kokkos::View<double*>& inv_x) const;
inline double acoustic_dt(const Kokkos::View<const double*>& Vj,
const Kokkos::View<const double*>& cj) const;
- inline void computeExplicitFluxes(const Kokkos::View<const double*>& xr,
- const Kokkos::View<const double*>& xj,
+ inline void computeExplicitFluxes(const Kokkos::View<const Rd*>& xr,
+ const Kokkos::View<const Rd*>& xj,
const Kokkos::View<const double*>& rhoj,
- const Kokkos::View<const double*>& uj,
+ const Kokkos::View<const Rd*>& uj,
const Kokkos::View<const double*>& pj,
const Kokkos::View<const double*>& cj,
const Kokkos::View<const double*>& Vj,
- const Kokkos::View<const double*[2]>& Cjr,
+ const Kokkos::View<const Rd*[2]>& Cjr,
const Kokkos::View<const int*[2]>& cell_nodes,
const Kokkos::View<const int*[2]>& node_cells,
const Kokkos::View<const int*>& node_nb_cells,
const Kokkos::View<const int*[2]>& node_cell_local_node,
- Kokkos::View<double*>& ur,
- Kokkos::View<double*[2]>& Fjr);
+ Kokkos::View<Rd*>& ur,
+ Kokkos::View<Rd*[2]>& Fjr);
struct ReduceMin;
const size_t m_nj;
const size_t m_nr;
- Kokkos::View<double*> m_br;
- Kokkos::View<double*[2]> m_Ajr;
- Kokkos::View<double*> m_Ar;
- Kokkos::View<double*> m_inv_Ar;
- Kokkos::View<double*[2]> m_Fjr;
- Kokkos::View<double*> m_ur;
+ Kokkos::View<Rd*> m_br;
+ Kokkos::View<Rdd*[2]> m_Ajr;
+ Kokkos::View<Rdd*> m_Ar;
+ Kokkos::View<Rdd*> m_inv_Ar;
+ Kokkos::View<Rd*[2]> m_Fjr;
+ Kokkos::View<Rd*> m_ur;
Kokkos::View<double*> m_rhocj;
Kokkos::View<double*> m_Vj_over_cj;
diff --git a/main.cpp b/main.cpp
index 8db2c2a5fe24f5c991d92a6fd17ca106498cafe2..185c41dbd5cafaa86b2d312819532c410ff28a83 100644
--- a/main.cpp
+++ b/main.cpp
@@ -12,6 +12,7 @@
#include <AcousticSolverTest.hpp>
#include <TinyVector.hpp>
+#include <TinyMatrix.hpp>
#include <CLI/CLI.hpp>
#include <cassert>
@@ -81,19 +82,19 @@ int main(int argc, char *argv[])
std::map<std::string, double> method_cost_map;
- { // Basic function based acoustic solver
- Kokkos::Timer timer;
- timer.reset();
- RawKokkos::AcousticSolver(number);
- method_cost_map["RawKokkos"] = timer.seconds();
- }
+ // { // Basic function based acoustic solver
+ // Kokkos::Timer timer;
+ // timer.reset();
+ // RawKokkos::AcousticSolver(number);
+ // method_cost_map["RawKokkos"] = timer.seconds();
+ // }
- { // class for acoustic solver (mesh less)
- Kokkos::Timer timer;
- timer.reset();
- MeshLessAcousticSolver acoustic_solver(number);
- method_cost_map["MeshLessAcousticSolver"] = timer.seconds();
- }
+ // { // class for acoustic solver (mesh less)
+ // Kokkos::Timer timer;
+ // timer.reset();
+ // MeshLessAcousticSolver acoustic_solver(number);
+ // method_cost_map["MeshLessAcousticSolver"] = timer.seconds();
+ // }
{ // class for acoustic solver test
Kokkos::Timer timer;
@@ -102,12 +103,12 @@ int main(int argc, char *argv[])
method_cost_map["AcousticSolverTest"] = timer.seconds();
}
- { // class for acoustic solver
- Kokkos::Timer timer;
- timer.reset();
- AcousticSolver acoustic_solver(number);
- method_cost_map["AcousticSolver"] = timer.seconds();
- }
+ // { // class for acoustic solver
+ // Kokkos::Timer timer;
+ // timer.reset();
+ // AcousticSolver acoustic_solver(number);
+ // method_cost_map["AcousticSolver"] = timer.seconds();
+ // }
Kokkos::View<TinyVector<2,double>*[2]> test("test", 10);
constexpr size_t N = 3;
@@ -136,6 +137,9 @@ int main(int argc, char *argv[])
TinyVector<2> x=zero;
TinyVector<2> y{1,2};
+
+ TinyVector<6> kk{1,2,3,4,5,6};
+
x=TinyVector<2>{3,2};
std::cout << x << "-" << y << "=" << x-y << std::endl;
@@ -143,5 +147,25 @@ int main(int argc, char *argv[])
std::cout << "3*" << x << '=' << 3*x << std::endl;
x[1]=3; y[1]=0.1;
std::cout << "< " << x << " | " << y << " > = " << (x,y) << std::endl;
+
+
+ TinyMatrix<2> Z = zero;
+ TinyMatrix<2> A(1,2,
+ 3,4);
+ TinyMatrix<2> B = identity;
+
+ B(1,1)=2;
+
+ std::cout << "A=" << A << '\n';
+ std::cout << "B=" << B << '\n';
+ std::cout << "2A=" << 2*A << '\n';
+ std::cout << "A+B=" << A+B << '\n';
+ std::cout << "A-B=" << A-B << '\n';
+ std::cout << "AB=" << A*B << '\n';
+
+ std::cout << "Ax=" << A*x << '\n';
+
+ std::cout << "x tens y" << tensorProduct(x,y) << '\n';
+
return 0;
}