diff --git a/src/mesh/Connectivity1D.hpp b/src/mesh/Connectivity1D.hpp index 21d89450c7ea14d0d9ac8f8467c5e929d98a057f..c6cb77aa17130d3e98c9e55a5792cfb361bc721c 100644 --- a/src/mesh/Connectivity1D.hpp +++ b/src/mesh/Connectivity1D.hpp @@ -95,21 +95,11 @@ public: return m_max_nb_node_per_cell; } - const Kokkos::View<const unsigned int**> cellNodes() const - { - return m_cell_nodes; - } - const Kokkos::View<const unsigned int**> cellFaces() const { return m_cell_faces; } - const Kokkos::View<const unsigned short*> cellNbNodes() const - { - return m_cell_nb_nodes; - } - const Kokkos::View<const double*> invCellNbNodes() const { return m_inv_cell_nb_nodes; diff --git a/src/mesh/Connectivity2D.hpp b/src/mesh/Connectivity2D.hpp index 0bab4bf58d349575eb42f2baa920bce85f71c382..0e5cb2685cd8218527aed0ad6f1c38b2ed8ac855 100644 --- a/src/mesh/Connectivity2D.hpp +++ b/src/mesh/Connectivity2D.hpp @@ -213,31 +213,11 @@ class Connectivity2D return m_max_nb_node_per_cell; } - const Kokkos::View<const unsigned int**> cellNodes() const - { - return m_cell_nodes; - } - - // const Kokkos::View<const unsigned int**> cellFaces() const - // { - // return m_cell_faces; - // } - - const Kokkos::View<const unsigned short*> cellNbNodes() const - { - return m_cell_nb_nodes; - } - const Kokkos::View<const double*> invCellNbNodes() const { return m_inv_cell_nb_nodes; } - // const Kokkos::View<const unsigned short*> cellNbFaces() const - // { - // return m_cell_nb_faces; - // } - const Kokkos::View<const unsigned short**> faceCellLocalFace() const { return m_face_cell_local_face; diff --git a/src/mesh/Connectivity3D.hpp b/src/mesh/Connectivity3D.hpp index 3a92bcfc9f3b88b17cdc6d550cf9677020380c1f..ee1b252ebf7f969d0994832b061ac00f8a074896 100644 --- a/src/mesh/Connectivity3D.hpp +++ b/src/mesh/Connectivity3D.hpp @@ -426,14 +426,16 @@ private: return m_node_to_cell_matrix.numRows(); } + KOKKOS_INLINE_FUNCTION const size_t& numberOfFaces() const { return m_number_of_faces; } - const size_t& numberOfCells() const + KOKKOS_INLINE_FUNCTION + size_t numberOfCells() const { - return m_number_of_cells; + return m_cell_to_node_matrix.numRows(); } const size_t& maxNbFacePerCell() const @@ -451,16 +453,6 @@ private: return m_max_nb_node_per_face; } - const Kokkos::View<const unsigned int**> cellNodes() const - { - return m_cell_nodes; - } - - const Kokkos::View<const unsigned short*> cellNbNodes() const - { - return m_cell_nb_nodes; - } - const Kokkos::View<const double*> invCellNbNodes() const { return m_inv_cell_nb_nodes; @@ -492,6 +484,7 @@ private: } m_cell_nb_nodes = cell_nb_nodes; } + { Kokkos::View<unsigned int**> cell_nodes("cell_nodes", cell_by_node_vector.size(), 8); for (size_t j=0; j<cell_by_node_vector.size(); ++j) { @@ -501,6 +494,7 @@ private: } m_cell_nodes = cell_nodes; } + { Kokkos::View<double*> inv_cell_nb_nodes("inv_cell_nb_nodes", m_number_of_cells); Kokkos::parallel_for(m_number_of_cells, KOKKOS_LAMBDA(const int& j){ diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp index 11516ed001bf8bf182ca55c4f00218f183770a31..f0a4da3c66e4824a552cac6995f7b96a4e65d603 100644 --- a/src/mesh/GmshReader.cpp +++ b/src/mesh/GmshReader.cpp @@ -879,29 +879,6 @@ GmshReader::__proceedData() } } - // const Kokkos::View<unsigned int**> cell_nodes("cell_nodes", nb_cells, max_nb_node_per_cell); - - // for (size_t j=0; j<nb_triangles; ++j) { - // cell_nodes(j,0) = __triangles[j][0]; - // cell_nodes(j,1) = __triangles[j][1]; - // cell_nodes(j,2) = __triangles[j][2]; - // } - - // const size_t nb_quadrangles = __quadrangles.extent(0); - // for (size_t j=0; j<nb_quadrangles; ++j) { - // const size_t jq = j+nb_triangles; - // cell_nodes(jq,0) = __quadrangles[j][0]; - // cell_nodes(jq,1) = __quadrangles[j][1]; - // cell_nodes(jq,2) = __quadrangles[j][2]; - // cell_nodes(jq,3) = __quadrangles[j][3]; - // } - // const Kokkos::View<unsigned short*> cell_nb_nodes("cell_nb_nodes", nb_cells); - // for (size_t j=0; j<nb_triangles; ++j) { - // cell_nb_nodes[j] = 3; - // } - // for (size_t j=nb_triangles; j<nb_triangles+nb_quadrangles; ++j) { - // cell_nb_nodes[j] = 4; - // } std::shared_ptr<Connectivity2D> p_connectivity(new Connectivity2D(cell_by_node_vector)); Connectivity2D& connectivity = *p_connectivity; diff --git a/src/mesh/Mesh.hpp b/src/mesh/Mesh.hpp index 9d3d596c26dd5117890e1071771f47fba0057d08..3aee2c287d5cdd06c4a754115caf930973c4f2cb 100644 --- a/src/mesh/Mesh.hpp +++ b/src/mesh/Mesh.hpp @@ -42,12 +42,14 @@ public: return m_connectivity->numberOfNodes(); } - const size_t& numberOfFaces() const + KOKKOS_INLINE_FUNCTION + size_t numberOfFaces() const { return m_connectivity->numberOfFaces(); } - const size_t& numberOfCells() const + KOKKOS_INLINE_FUNCTION + size_t numberOfCells() const { return m_connectivity->numberOfCells(); } diff --git a/src/mesh/MeshData.hpp b/src/mesh/MeshData.hpp index e2c0e0ac612cf46161170300dc4161c19accb8ab..f9732fbe98cec1c312b2d244a33e964493f1a9be 100644 --- a/src/mesh/MeshData.hpp +++ b/src/mesh/MeshData.hpp @@ -30,45 +30,39 @@ private: KOKKOS_INLINE_FUNCTION void _updateCenter() { // Computes vertices isobarycenter - if(dimension == 1) { + if constexpr (dimension == 1) { const Kokkos::View<const Rd*> xr = m_mesh.xr(); - const Kokkos::View<const unsigned int**>& cell_nodes - = m_mesh.connectivity().cellNodes(); + Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ - m_xj[j] = 0.5*(xr[cell_nodes(j,0)]+xr[cell_nodes(j,1)]); - }); + const auto& cell_nodes = m_mesh.connectivity().m_cell_to_node_matrix.rowConst(j); + m_xj[j] = 0.5*(xr[cell_nodes(0)]+xr[cell_nodes(1)]); + }); } else { const Kokkos::View<const Rd*> xr = m_mesh.xr(); - const Kokkos::View<const unsigned int**>& cell_nodes - = m_mesh.connectivity().cellNodes(); - const Kokkos::View<const unsigned short*>& cell_nb_nodes - = m_mesh.connectivity().cellNbNodes(); const Kokkos::View<const double*>& inv_cell_nb_nodes - = m_mesh.connectivity().invCellNbNodes(); + = m_mesh.connectivity().invCellNbNodes(); Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ - Rd X = zero; - for (int R=0; R<cell_nb_nodes[j]; ++R) { - X += xr[cell_nodes(j,R)]; - } - m_xj[j] = inv_cell_nb_nodes[j]*X; - }); + Rd X = zero; + const auto& cell_nodes = m_mesh.connectivity().m_cell_to_node_matrix.rowConst(j); + for (size_t R=0; R<cell_nodes.length; ++R) { + X += xr[cell_nodes(R)]; + } + m_xj[j] = inv_cell_nb_nodes[j]*X; + }); } } KOKKOS_INLINE_FUNCTION void _updateVolume() { - const Kokkos::View<const unsigned int**>& cell_nodes - = m_mesh.connectivity().cellNodes(); - const Kokkos::View<const unsigned short*> cell_nb_nodes - = m_mesh.connectivity().cellNbNodes(); - const Kokkos::View<const Rd*> xr = m_mesh.xr(); Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ double sum_cjr_xr = 0; - for (int R=0; R<cell_nb_nodes[j]; ++R) { - sum_cjr_xr += (xr[cell_nodes(j,R)], m_Cjr(j,R)); + const auto& cell_nodes = m_mesh.connectivity().m_cell_to_node_matrix.rowConst(j); + + for (size_t R=0; R<cell_nodes.length; ++R) { + sum_cjr_xr += (xr[cell_nodes(R)], m_Cjr(j,R)); } m_Vj[j] = inv_dimension * sum_cjr_xr; }); @@ -80,42 +74,35 @@ private: // Cjr/njr/ljr are constant overtime } else if constexpr (dimension == 2) { - const Kokkos::View<const unsigned int**>& cell_nodes - = m_mesh.connectivity().cellNodes(); - const Kokkos::View<const unsigned short*> cell_nb_nodes - = m_mesh.connectivity().cellNbNodes(); - const Kokkos::View<const Rd*> xr = m_mesh.xr(); Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ - for (int R=0; R<cell_nb_nodes[j]; ++R) { - int Rp1 = (R+1)%cell_nb_nodes[j]; - int Rm1 = (R+cell_nb_nodes[j]-1)%cell_nb_nodes[j]; - Rd half_xrp_xrm = 0.5*(xr(cell_nodes(j,Rp1))-xr(cell_nodes(j,Rm1))); + const auto& cell_nodes = m_mesh.connectivity().m_cell_to_node_matrix.rowConst(j); + for (size_t R=0; R<cell_nodes.length; ++R) { + int Rp1 = (R+1)%cell_nodes.length; + int Rm1 = (R+cell_nodes.length-1)%cell_nodes.length; + Rd half_xrp_xrm = 0.5*(xr(cell_nodes(Rp1))-xr(cell_nodes(Rm1))); m_Cjr(j,R) = Rd{-half_xrp_xrm[1], half_xrp_xrm[0]}; } }); Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ - for (int R=0; R<cell_nb_nodes[j]; ++R) { + const auto& cell_nodes = m_mesh.connectivity().m_cell_to_node_matrix.rowConst(j); + for (size_t R=0; R<cell_nodes.length; ++R) { const Rd& Cjr = m_Cjr(j,R); m_ljr(j,R) = l2Norm(Cjr); } }); Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ - for (int R=0; R<cell_nb_nodes[j]; ++R) { + const auto& cell_nodes = m_mesh.connectivity().m_cell_to_node_matrix.rowConst(j); + for (size_t R=0; R<cell_nodes.length; ++R) { const Rd& Cjr = m_Cjr(j,R); const double inv_ljr = 1./m_ljr(j,R); m_njr(j,R) = inv_ljr*Cjr; } }); } else if (dimension ==3) { - const Kokkos::View<const unsigned int**>& cell_nodes - = m_mesh.connectivity().cellNodes(); - const Kokkos::View<const unsigned short*> cell_nb_nodes - = m_mesh.connectivity().cellNbNodes(); - const Kokkos::View<const Rd*> xr = m_mesh.xr(); Kokkos::View<Rd**> Nlr("Nlr", m_mesh.connectivity().numberOfFaces(), m_mesh.connectivity().maxNbNodePerFace()); @@ -141,12 +128,13 @@ private: }); Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) { - for (int R=0; R<cell_nb_nodes[j]; ++R) { + const auto& cell_nodes = m_mesh.connectivity().m_cell_to_node_matrix.rowConst(j); + for (size_t R=0; R<cell_nodes.length; ++R) { m_Cjr(j,R) = zero; } std::map<unsigned int, unsigned short> node_id_to_local; - for (size_t R=0; R<cell_nb_nodes[j]; ++R) { - node_id_to_local[cell_nodes(j,R)] = R; + for (size_t R=0; R<cell_nodes.length; ++R) { + node_id_to_local[cell_nodes(R)] = R; } const auto& cell_faces = m_mesh.connectivity().m_cell_to_face_matrix.rowConst(j); const auto& cell_faces_is_reversed = m_mesh.connectivity().m_cell_to_face_is_reversed_matrix.rowConst(j); @@ -166,14 +154,16 @@ private: }); Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) { - for (int R=0; R<cell_nb_nodes[j]; ++R) { + const auto& cell_nodes = m_mesh.connectivity().m_cell_to_node_matrix.rowConst(j); + for (size_t R=0; R<cell_nodes.length; ++R) { const Rd& Cjr = m_Cjr(j,R); m_ljr(j,R) = l2Norm(Cjr); } }); Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) { - for (int R=0; R<cell_nb_nodes[j]; ++R) { + const auto& cell_nodes = m_mesh.connectivity().m_cell_to_node_matrix.rowConst(j); + for (size_t R=0; R<cell_nodes.length; ++R) { const Rd& Cjr = m_Cjr(j,R); const double inv_ljr = 1./m_ljr(j,R); m_njr(j,R) = inv_ljr*Cjr; diff --git a/src/mesh/MeshNodeBoundary.hpp b/src/mesh/MeshNodeBoundary.hpp index 4349a0a5b67dc372030603f17dd3d680dde6eae4..9785552fe22d0f937bf4e05f6fc8c197ec8452d7 100644 --- a/src/mesh/MeshNodeBoundary.hpp +++ b/src/mesh/MeshNodeBoundary.hpp @@ -305,14 +305,13 @@ _getOutgoingNormal(const MeshType& mesh) const R normal = this->_getNormal(mesh); const Kokkos::View<const R*>& xr = mesh.xr(); - const Kokkos::View<const unsigned int**>& cell_nodes = mesh.connectivity().cellNodes(); - const Kokkos::View<const unsigned short*>& cell_nb_nodes = mesh.connectivity().cellNbNodes(); const size_t r0 = m_node_list[0]; const size_t j0 = mesh.connectivity().m_node_to_cell_matrix.rowConst(r0)(0); + const auto& j0_nodes = mesh.connectivity().m_cell_to_node_matrix.rowConst(j0); double max_height = 0; - for (int r=0; r<cell_nb_nodes(j0); ++r) { - const double height = (xr(cell_nodes(j0, r))-xr(r0), normal); + for (size_t r=0; r<j0_nodes.length; ++r) { + const double height = (xr(j0_nodes(r))-xr(r0), normal); if (std::abs(height) > std::abs(max_height)) { max_height = height; } @@ -336,14 +335,13 @@ _getOutgoingNormal(const MeshType& mesh) const R2 normal = this->_getNormal(mesh); const Kokkos::View<const R2*>& xr = mesh.xr(); - const Kokkos::View<const unsigned int**>& cell_nodes = mesh.connectivity().cellNodes(); - const Kokkos::View<const unsigned short*>& cell_nb_nodes = mesh.connectivity().cellNbNodes(); const size_t r0 = m_node_list[0]; const size_t j0 = mesh.connectivity().m_node_to_cell_matrix.rowConst(r0)(0); + const auto& j0_nodes = mesh.connectivity().m_cell_to_node_matrix.rowConst(j0); double max_height = 0; - for (int r=0; r<cell_nb_nodes(j0); ++r) { - const double height = (xr(cell_nodes(j0, r))-xr(r0), normal); + for (size_t r=0; r<j0_nodes.length; ++r) { + const double height = (xr(j0_nodes(r))-xr(r0), normal); if (std::abs(height) > std::abs(max_height)) { max_height = height; } @@ -367,14 +365,13 @@ _getOutgoingNormal(const MeshType& mesh) const R3 normal = this->_getNormal(mesh); const Kokkos::View<const R3*>& xr = mesh.xr(); - const Kokkos::View<const unsigned int**>& cell_nodes = mesh.connectivity().cellNodes(); - const Kokkos::View<const unsigned short*>& cell_nb_nodes = mesh.connectivity().cellNbNodes(); const size_t r0 = m_node_list[0]; const size_t j0 = mesh.connectivity().m_node_to_cell_matrix.rowConst(r0)(0); + const auto& j0_nodes = mesh.connectivity().m_cell_to_node_matrix.rowConst(j0); double max_height = 0; - for (int r=0; r<cell_nb_nodes(j0); ++r) { - const double height = (xr(cell_nodes(j0, r))-xr(r0), normal); + for (size_t r=0; r<j0_nodes.length; ++r) { + const double height = (xr(j0_nodes(r))-xr(r0), normal); if (std::abs(height) > std::abs(max_height)) { max_height = height; } diff --git a/src/output/VTKWriter.hpp b/src/output/VTKWriter.hpp index 4e34f666cc0829eef509f4140aaa3aab6d9f9727..f297325bd328a15cb57c144021dcf5d9b5c7e302 100644 --- a/src/output/VTKWriter.hpp +++ b/src/output/VTKWriter.hpp @@ -67,13 +67,12 @@ class VTKWriter fout << "</Points>\n"; fout << "<Cells>\n"; - const Kokkos::View<const unsigned int**> cell_nodes = mesh.connectivity().cellNodes(); - const Kokkos::View<const unsigned short*> cell_nb_nodes = mesh.connectivity().cellNbNodes(); fout << "<DataArray type=\"Int32\" Name=\"connectivity\" NumberOfComponents=\"1\" format=\"ascii\">\n"; for (unsigned int j=0; j<mesh.numberOfCells(); ++j) { - for (unsigned short r=0; r<cell_nb_nodes[j]; ++r) { - fout << cell_nodes(j,r) << ' '; + const auto& cell_nodes = mesh.connectivity().m_cell_to_node_matrix.rowConst(j); + for (unsigned short r=0; r<cell_nodes.length; ++r) { + fout << cell_nodes(r) << ' '; } } fout << '\n'; @@ -83,7 +82,8 @@ class VTKWriter { unsigned int offset=0; for (unsigned int j=0; j<mesh.numberOfCells(); ++j) { - offset += cell_nb_nodes[j]; + const auto& cell_nodes = mesh.connectivity().m_cell_to_node_matrix.rowConst(j); + offset += cell_nodes.length; fout << offset << ' '; } } @@ -92,7 +92,8 @@ class VTKWriter fout << "<DataArray type=\"Int8\" Name=\"types\" NumberOfComponents=\"1\" format=\"ascii\">\n"; for (unsigned int j=0; j<mesh.numberOfCells(); ++j) { - switch (cell_nb_nodes[j]) { + const auto& cell_nodes = mesh.connectivity().m_cell_to_node_matrix.rowConst(j); + switch (cell_nodes.length) { case 2: { fout << "3 "; break; diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp index dadcc82fbd31b916455593e9e099f2218da52a60..0ce9beef810467a35466247c04a1abccb233aa1e 100644 --- a/src/scheme/AcousticSolver.hpp +++ b/src/scheme/AcousticSolver.hpp @@ -82,12 +82,12 @@ private: const Kokkos::View<const Rdd**> computeAjr(const Kokkos::View<const double*>& rhocj, const Kokkos::View<const double**>& ljr, - const Kokkos::View<const Rd**>& njr) { - const Kokkos::View<const unsigned short*> cell_nb_nodes - = m_connectivity.cellNbNodes(); - + const Kokkos::View<const Rd**>& njr) + { Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) { - for (int r=0; r<cell_nb_nodes[j]; ++r) { + const auto& cell_nodes = m_mesh.connectivity().m_cell_to_node_matrix.rowConst(j); + + for (size_t r=0; r<cell_nodes.length; ++r) { m_Ajr(j,r) = tensorProduct((rhocj(j)*ljr(j,r))*njr(j,r), njr(j,r)); } }); @@ -179,7 +179,8 @@ private: Kokkos::View<Rd*> computeUr(const Kokkos::View<const Rdd*>& Ar, - const Kokkos::View<const Rd*>& br) { + const Kokkos::View<const Rd*>& br) + { inverse(Ar, m_inv_Ar); const Kokkos::View<const Rdd*> invAr = m_inv_Ar; Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const int& r) { @@ -194,14 +195,13 @@ private: const Kokkos::View<const Rd*>& ur, const Kokkos::View<const Rd**>& Cjr, const Kokkos::View<const Rd*>& uj, - const Kokkos::View<const double*>& pj) { - const Kokkos::View<const unsigned int**>& cell_nodes = m_connectivity.cellNodes(); - const Kokkos::View<const unsigned short*> cell_nb_nodes - = m_connectivity.cellNbNodes(); - + const Kokkos::View<const double*>& pj) + { Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) { - for (int r=0; r<cell_nb_nodes[j]; ++r) { - m_Fjr(j,r) = Ajr(j,r)*(uj(j)-ur(cell_nodes(j,r)))+pj(j)*Cjr(j,r); + const auto& cell_nodes = m_mesh.connectivity().m_cell_to_node_matrix.rowConst(j); + + for (size_t r=0; r<cell_nodes.length; ++r) { + m_Fjr(j,r) = Ajr(j,r)*(uj(j)-ur(cell_nodes(r)))+pj(j)*Cjr(j,r); } }); @@ -209,7 +209,8 @@ private: } void inverse(const Kokkos::View<const Rdd*>& A, - Kokkos::View<Rdd*>& inv_A) const { + Kokkos::View<Rdd*>& inv_A) const + { Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const int& r) { inv_A(r) = ::inverse(A(r)); }); @@ -271,14 +272,15 @@ public: KOKKOS_INLINE_FUNCTION double acoustic_dt(const Kokkos::View<const double*>& Vj, - const Kokkos::View<const double*>& cj) const { + const Kokkos::View<const double*>& cj) const + { const Kokkos::View<const double**> ljr = m_mesh_data.ljr(); - const Kokkos::View<const unsigned short*>& cell_nb_nodes - = m_connectivity.cellNbNodes(); Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ + const auto& cell_nodes = m_mesh.connectivity().m_cell_to_node_matrix.rowConst(j); + double S = 0; - for (int r=0; r<cell_nb_nodes(j); ++r) { + for (size_t r=0; r<cell_nodes.length; ++r) { S += ljr(j,r); } m_Vj_over_cj[j] = 2*Vj[j]/(S*cj[j]); @@ -314,17 +316,15 @@ public: const Kokkos::View<const Rd**> Fjr = m_Fjr; const Kokkos::View<const Rd*> ur = m_ur; - const Kokkos::View<const unsigned short*> cell_nb_nodes - = m_connectivity.cellNbNodes(); - const Kokkos::View<const unsigned int**>& cell_nodes - = m_connectivity.cellNodes(); const Kokkos::View<const double*> inv_mj = unknowns.invMj(); Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) { + const auto& cell_nodes = m_mesh.connectivity().m_cell_to_node_matrix.rowConst(j); + Rd momentum_fluxes = zero; double energy_fluxes = 0; - for (int R=0; R<cell_nb_nodes[j]; ++R) { - const int r=cell_nodes(j,R); + for (size_t R=0; R<cell_nodes.length; ++R) { + const unsigned int r=cell_nodes(R); momentum_fluxes += Fjr(j,R); energy_fluxes += (Fjr(j,R), ur[r]); }