diff --git a/.github/workflows/idefix-ci-jobs.yml b/.github/workflows/idefix-ci-jobs.yml index c0e62296..dfb89e0f 100644 --- a/.github/workflows/idefix-ci-jobs.yml +++ b/.github/workflows/idefix-ci-jobs.yml @@ -201,3 +201,5 @@ jobs: run: scripts/ci/run-tests $IDEFIX_DIR/test/utils/lookupTable -all $TESTME_OPTIONS - name: Dump Image run: scripts/ci/run-tests $IDEFIX_DIR/test/utils/dumpImage -all $TESTME_OPTIONS + - name: Column density + run: scripts/ci/run-tests $IDEFIX_DIR/test/utils/columnDensity -all $TESTME_OPTIONS diff --git a/doc/source/faq.rst b/doc/source/faq.rst index 75a07a68..4c7bc056 100644 --- a/doc/source/faq.rst +++ b/doc/source/faq.rst @@ -73,3 +73,6 @@ I want to test a modification to *Idefix* source code specific to my problem wit I want to use a lookup table from a CSV file in my idefix_loop. How could I proceed? Use the ``LookupTable`` class (see :ref:`LookupTableClass`) + +I want to compute a cumulative sum (e.g a column density) on the fly. How could I proceed? + Use the ``Column`` class (see :class:`::Column`) diff --git a/src/utils/CMakeLists.txt b/src/utils/CMakeLists.txt index b36f8a08..afb4f57d 100644 --- a/src/utils/CMakeLists.txt +++ b/src/utils/CMakeLists.txt @@ -5,4 +5,6 @@ target_sources(idefix PUBLIC ${CMAKE_CURRENT_LIST_DIR}/dumpImage.cpp PUBLIC ${CMAKE_CURRENT_LIST_DIR}/dumpImage.hpp PUBLIC ${CMAKE_CURRENT_LIST_DIR}/lookupTable.hpp + PUBLIC ${CMAKE_CURRENT_LIST_DIR}/column.cpp + PUBLIC ${CMAKE_CURRENT_LIST_DIR}/column.hpp ) diff --git a/src/utils/column.cpp b/src/utils/column.cpp new file mode 100644 index 00000000..34d55f56 --- /dev/null +++ b/src/utils/column.cpp @@ -0,0 +1,240 @@ +// *********************************************************************************** +// Idefix MHD astrophysical code +// Copyright(C) Geoffroy R. J. Lesur +// and other code contributors +// Licensed under CeCILL 2.1 License, see COPYING for more information +// *********************************************************************************** + + +#include +#include + +#include "column.hpp" +#include "idefix.hpp" +#include "input.hpp" +#include "output.hpp" +#include "grid.hpp" +#include "dataBlock.hpp" +#include "dataBlockHost.hpp" + +Column::Column(int dir, int sign, DataBlock *data) + : direction(dir), sign(sign) { + idfx::pushRegion("Column::Column"); + this->np_tot = data->np_tot; + this->np_int = data->np_int; + this->beg = data->beg; + + if(dir>= DIMENSIONS || dir < IDIR) IDEFIX_ERROR("Unknown direction for Column constructor"); + + // Allocate the array on which we do the average + this->ColumnArray = IdefixArray3D("ColumnArray",np_tot[KDIR], np_tot[JDIR], np_tot[IDIR]); + + // Elementary volumes and area + this->Volume = data->dV; + this->Area = data->A[dir]; + + // allocate helper array + if(dir == IDIR) { + localSum = IdefixArray2D("localSum",np_tot[KDIR], np_tot[JDIR]); + } + if(dir == JDIR) { + localSum = IdefixArray2D("localSum",np_tot[KDIR], np_tot[IDIR]); + } + if(dir == KDIR) { + localSum = IdefixArray2D("localSum",np_tot[JDIR], np_tot[IDIR]); + } + #ifdef WITH_MPI + // Create sub-MPI communicator dedicated to scan + int remainDims[3] = {false, false, false}; + remainDims[dir] = true; + MPI_Cart_sub(data->mygrid->CartComm, remainDims, &this->ColumnComm); + MPI_Comm_rank(this->ColumnComm, &this->MPIrank); + MPI_Comm_size(this->ColumnComm, &this->MPIsize); + + // create MPI class for boundary Xchanges + int ntarget = 0; + std::vector mapVars; + mapVars.push_back(ntarget); + + this->mpi.Init(data->mygrid, mapVars, data->nghost.data(), data->np_int.data()); + this->nproc = data->mygrid->nproc; + #endif + idfx::popRegion(); +} + +void Column::ComputeColumn(IdefixArray4D in, const int var) { + idfx::pushRegion("Column::ComputeColumn"); + const int nk = np_int[KDIR]; + const int nj = np_int[JDIR]; + const int ni = np_int[IDIR]; + const int kb = beg[KDIR]; + const int jb = beg[JDIR]; + const int ib = beg[IDIR]; + const int ke = kb+nk; + const int je = jb+nj; + const int ie = ib+ni; + + const int direction = this->direction; + auto column = this->ColumnArray; + auto dV = this->Volume; + auto A = this->Area; + auto localSum = this->localSum; + + if(direction==IDIR) { + // Inspired from loop.hpp + Kokkos::parallel_for("ColumnX1", team_policy (nk*nj, Kokkos::AUTO), + KOKKOS_LAMBDA (member_type team_member) { + int k = team_member.league_rank() / nj; + int j = team_member.league_rank() - k*nj + jb; + k += kb; + Kokkos::parallel_scan(Kokkos::TeamThreadRange<>(team_member,ib,ie), + [=] (int i, real &partial_sum, bool is_final) { + partial_sum += in(var,k,j,i)*dV(k,j,i) / (0.5*(A(k,j,i)+A(k,j,i+1))); + if(is_final) column(k,j,i) = partial_sum; + //partial_sum += in(var,k,j,i)*dV(k,j,i) / (0.5*(A(k,j,i)+A(k,j,i+1))); + }); + }); + } + if(direction==JDIR) { + // Inspired from loop.hpp + Kokkos::parallel_for("ColumnX2", team_policy (nk*ni, Kokkos::AUTO), + KOKKOS_LAMBDA (member_type team_member) { + int k = team_member.league_rank() / ni; + int i = team_member.league_rank() - k*ni + ib; + k += kb; + Kokkos::parallel_scan(Kokkos::TeamThreadRange<>(team_member,jb,je), + [=] (int j, real &partial_sum, bool is_final) { + partial_sum += in(var,k,j,i)*dV(k,j,i) / (0.5*(A(k,j,i)+A(k,j+1,i))); + if(is_final) column(k,j,i) = partial_sum; + }); + }); + } + if(direction==KDIR) { + // Inspired from loop.hpp + Kokkos::parallel_for("ColumnX2", team_policy (nj*ni, Kokkos::AUTO), + KOKKOS_LAMBDA (member_type team_member) { + int j = team_member.league_rank() / ni; + int i = team_member.league_rank() - j*ni + ib; + j += jb; + Kokkos::parallel_scan(Kokkos::TeamThreadRange<>(team_member,kb,ke), + [=] (int k, real &partial_sum, bool is_final) { + partial_sum += in(var,k,j,i)*dV(k,j,i) / (0.5*(A(k,j,i)+A(k+1,j,i))); + if(is_final) column(k,j,i) = partial_sum; + }); + }); + } + + #ifdef WITH_MPI + // Load the current sum + int dst,src; + MPI_Cart_shift(this->ColumnComm,0, 1, &src, &dst); + int size = localSum.extent(0)*localSum.extent(1); + if(MPIrank>0) { + MPI_Status status; + // Get the cumulative sum from previous processes + Kokkos::fence(); + MPI_Recv(localSum.data(), size, realMPI, src, 20, ColumnComm, &status); + // Add this to our cumulative sum + idefix_for("Addsum",kb,ke,jb,je,ib,ie, + KOKKOS_LAMBDA(int k, int j, int i) { + if(direction == IDIR) column(k,j,i) += localSum(k,j); + if(direction == JDIR) column(k,j,i) += localSum(k,i); + if(direction == KDIR) column(k,j,i) += localSum(j,i); + }); + } // rank =0 + // Get the final sum + if(MPIrank < MPIsize-1) { + if(direction==IDIR) { + idefix_for("Loadsum",kb,ke,jb,je, + KOKKOS_LAMBDA(int k, int j) { + localSum(k,j) = column(k,j,ie-1); + }); + } + if(direction==JDIR) { + idefix_for("Loadsum",kb,ke,ib,ie, + KOKKOS_LAMBDA(int k, int i) { + localSum(k,i) = column(k,je-1,i); + }); + } + if(direction==KDIR) { + idefix_for("Loadsum",jb,je,ib,ie, + KOKKOS_LAMBDA(int j, int i) { + localSum(j,i) = column(ke-1,j,i); + }); + } + // And send it + Kokkos::fence(); + MPI_Send(localSum.data(),size, realMPI, dst, 20, ColumnComm); + } // MPIrank small enough + #endif + // If we need it backwards + if(sign<0) { + // Compute total cumulative sum in the last subdomain + if(direction == IDIR) { + idefix_for("Loadsum",kb,ke,jb,je, + KOKKOS_LAMBDA(int k, int j) { + localSum(k,j) = column(k,j,ie-1) + in(var,k,j,ie-1) * dV(k,j,ie-1) + / (0.5*(A(k,j,ie-1)+A(k,j,ie))); + }); + } + if(direction == JDIR) { + idefix_for("Loadsum",kb,ke,ib,ie, + KOKKOS_LAMBDA(int k, int i) { + localSum(k,i) = column(k,je-1,i) + in(var,k,je-1,i) * dV(k,je-1,i) + / (0.5*(A(k,je-1,i)+A(k,je,i))); + }); + } + if(direction == KDIR) { + idefix_for("Loadsum",jb,je,ib,ie, + KOKKOS_LAMBDA(int j, int i) { + localSum(j,i) = column(ke-1,j,i) + in(var,ke-1,j,i) * dV(ke-1,j,i) + / (0.5*(A(ke-1,j,i)+A(ke,j,i))); + }); + } + + #ifdef WITH_MPI + Kokkos::fence(); + MPI_Bcast(localSum.data(),size, realMPI, MPIsize-1, ColumnComm); + #endif + // All substract the local column from the full column + + idefix_for("Subsum",kb,ke,jb,je,ib,ie, + KOKKOS_LAMBDA(int k, int j, int i) { + if(direction==IDIR) column(k,j,i) = localSum(k,j)-column(k,j,i); + if(direction==JDIR) column(k,j,i) = localSum(k,i)-column(k,j,i); + if(direction==KDIR) column(k,j,i) = localSum(j,i)-column(k,j,i); + }); + } // sign <0 + // Xchange boundary elements when using MPI to ensure that column + // density in the ghost zones are coherent + #ifdef WITH_MPI + // Create a 4D array that contains our column data + IdefixArray4D arr4D(column.data(), 1, this->np_tot[KDIR], + this->np_tot[JDIR], + this->np_tot[IDIR]); + + for(int dir = 0 ; dir < DIMENSIONS ; dir++) { + // MPI Exchange data when needed + if(nproc[dir]>1) { + switch(dir) { + case 0: + this->mpi.ExchangeX1(arr4D); + break; + case 1: + this->mpi.ExchangeX2(arr4D); + break; + case 2: + this->mpi.ExchangeX3(arr4D); + break; + } + } + } + #endif + idfx::popRegion(); +} + +void Column::ComputeColumn(IdefixArray3D in) { + // 4D alias + IdefixArray4D arr4D(in.data(), 1, in.extent(0), in.extent(1), in.extent(2)); + return this->ComputeColumn(arr4D,0); +} diff --git a/src/utils/column.hpp b/src/utils/column.hpp new file mode 100644 index 00000000..43383cbe --- /dev/null +++ b/src/utils/column.hpp @@ -0,0 +1,74 @@ +// *********************************************************************************** +// Idefix MHD astrophysical code +// Copyright(C) Geoffroy R. J. Lesur +// and other code contributors +// Licensed under CeCILL 2.1 License, see COPYING for more information +// *********************************************************************************** + +#ifndef UTILS_COLUMN_HPP_ +#define UTILS_COLUMN_HPP_ + +#include +#include + +#include "idefix.hpp" +#include "dataBlock.hpp" + + + +// A class to implement a parralel cumulative sum +class Column { + public: + //////////////////////////////////////////////////////////////////////////////////// + /// @brief Constructor of a cumulative sum(=column density) from setup's arguments + /// @param dir direction along which the integration is performed + /// @param sign: +1 for an integration from left to right, -1 for an integration from right to + /// left i.e (backwards) + /////////////////////////////////////////////////////////////////////////////////// + Column(int dir, int sign, DataBlock *); + + /////////////////////////////////////////////////////////////////////////////////// + /// @brief Effectively compute integral from the input array in argument + /// @param in: 4D input array + /// @param variable: index of the variable along which we do the integral (since the + /// intput array is 4D) + /////////////////////////////////////////////////////////////////////////////////// + void ComputeColumn(IdefixArray4D in, int variable); + + /////////////////////////////////////////////////////////////////////////////////// + /// @brief Effectively compute integral from the input array in argument + /// @param in: 3D input array + /////////////////////////////////////////////////////////////////////////////////// + void ComputeColumn(IdefixArray3D in); + + /////////////////////////////////////////////////////////////////////////////////// + /// @brief Get a reference to the computed column density array + /////////////////////////////////////////////////////////////////////////////////// + IdefixArray3D GetColumn() { + return (this->ColumnArray); + } + + private: + IdefixArray3D ColumnArray; + int direction; // The direction along which the column is computed + int sign; // whether we integrate from the left or from the right + std::array np_tot; + std::array np_int; + std::array beg; + + IdefixArray3D Area; + IdefixArray3D Volume; + + IdefixArray2D localSum; + #ifdef WITH_MPI + Mpi mpi; // Mpi object when WITH_MPI is set + MPI_Comm ColumnComm; + int MPIrank; + int MPIsize; + + std::array nproc; // 3D size of the MPI cartesian geometry + + #endif +}; + +#endif // UTILS_COLUMN_HPP_ diff --git a/test/utils/columnDensity/definitions.hpp b/test/utils/columnDensity/definitions.hpp new file mode 100644 index 00000000..3b581bd6 --- /dev/null +++ b/test/utils/columnDensity/definitions.hpp @@ -0,0 +1,4 @@ +#define COMPONENTS 3 +#define DIMENSIONS 3 + +#define GEOMETRY CARTESIAN diff --git a/test/utils/columnDensity/idefix.ini b/test/utils/columnDensity/idefix.ini new file mode 100644 index 00000000..5e80c8bc --- /dev/null +++ b/test/utils/columnDensity/idefix.ini @@ -0,0 +1,24 @@ +[Grid] +X1-grid 1 0.0 64 u 1.0 +X2-grid 1 0.0 32 u 1.0 +X3-grid 1 0.0 16 u 1.0 + +[TimeIntegrator] +CFL 0.8 +tstop 0.0 +nstages 2 + +[Hydro] +solver hllc +gamma 1.4 + +[Boundary] +X1-beg periodic +X1-end periodic +X2-beg periodic +X2-end periodic +X3-beg periodic +X3-end periodic + +[Output] +analysis 0.01 diff --git a/test/utils/columnDensity/setup.cpp b/test/utils/columnDensity/setup.cpp new file mode 100644 index 00000000..fe7c57e6 --- /dev/null +++ b/test/utils/columnDensity/setup.cpp @@ -0,0 +1,169 @@ +#include "idefix.hpp" +#include "setup.hpp" +#include "column.hpp" + + +Column *columnX1Left; +Column *columnX1Right; +Column *columnX2Left; +Column *columnX2Right; +Column *columnX3Left; +Column *columnX3Right; + +// Analyse data to check that column density works as expected +void Analysis(DataBlock & data) { + + DataBlockHost d(data); + + // Try the 4D array interface + columnX1Left->ComputeColumn(data.hydro->Vc,RHO); + columnX1Right->ComputeColumn(data.hydro->Vc,RHO); + + // Try the 3D array interface + IdefixArray3D rho("rho",data.np_tot[KDIR],data.np_tot[JDIR],data.np_tot[IDIR]); + auto Vc = data.hydro->Vc; + + idefix_for("init rho",data.beg[KDIR],data.end[KDIR],data.beg[JDIR],data.end[JDIR],data.beg[IDIR],data.end[IDIR], + KOKKOS_LAMBDA(int k, int j, int i) { + rho(k,j,i) = Vc(RHO,k,j,i); + }); + columnX2Left->ComputeColumn(rho); + columnX2Right->ComputeColumn(rho); + columnX3Left->ComputeColumn(rho); + columnX3Right->ComputeColumn(rho); + + IdefixArray3D columnDensityLeft, columnDensityRight; + IdefixArray3D::HostMirror columnDensityLeftHost, columnDensityRightHost; + // IDIR + columnDensityLeft = columnX1Left->GetColumn(); + columnDensityRight = columnX1Right->GetColumn(); + columnDensityLeftHost = Kokkos::create_mirror_view(columnDensityLeft); + columnDensityRightHost = Kokkos::create_mirror_view(columnDensityRight); + Kokkos::deep_copy(columnDensityLeftHost,columnDensityLeft); + Kokkos::deep_copy(columnDensityRightHost,columnDensityRight); + + real errMax = 0.0; + + // Because this particular setup has a constant density =1, we expect + // the column density from the left to match the cell right coordinates, + // and the column density from the right to match L-the cell left coordinate (where L is the box size) + // here we have L=1 + // This routine checks that all of the available column densities do match the expected values. + for(int k = data.beg[KDIR]; k < data.end[KDIR] ; k++) { + for(int j = data.beg[JDIR]; j < data.end[JDIR] ; j++) { + for(int i = data.beg[IDIR]; i < data.end[IDIR] ; i++) { + real err = std::fabs(columnDensityLeftHost(k,j,i)-d.xr[IDIR](i)); + if(err>errMax) errMax=err; + err = std::fabs(columnDensityRightHost(k,j,i)-(1-d.xl[IDIR](i))); + if(err>errMax) errMax=err; + } + } + } + idfx::cout << "Error on column density in IDIR=" << std::scientific << errMax << std::endl; + if(errMax>1e-14) { + IDEFIX_ERROR("Error above tolerance"); + } + // JDIR + columnDensityLeft = columnX2Left->GetColumn(); + columnDensityRight = columnX2Right->GetColumn(); + columnDensityLeftHost = Kokkos::create_mirror_view(columnDensityLeft); + columnDensityRightHost = Kokkos::create_mirror_view(columnDensityRight); + Kokkos::deep_copy(columnDensityLeftHost,columnDensityLeft); + Kokkos::deep_copy(columnDensityRightHost,columnDensityRight); + + errMax = 0.0; + + for(int k = data.beg[KDIR]; k < data.end[KDIR] ; k++) { + for(int j = data.beg[JDIR]; j < data.end[JDIR] ; j++) { + for(int i = data.beg[IDIR]; i < data.end[IDIR] ; i++) { + real err = std::fabs(columnDensityLeftHost(k,j,i)-d.xr[JDIR](j)); + if(err>errMax) errMax=err; + err = std::fabs(columnDensityRightHost(k,j,i)-(1-d.xl[JDIR](j))); + if(err>errMax) errMax=err; + } + } + } + idfx::cout << "Error on column density in JDIR=" << std::scientific << errMax << std::endl; + if(errMax>1e-14) { + IDEFIX_ERROR("Error above tolerance"); + } + // KDIR + columnDensityLeft = columnX3Left->GetColumn(); + columnDensityRight = columnX3Right->GetColumn(); + columnDensityLeftHost = Kokkos::create_mirror_view(columnDensityLeft); + columnDensityRightHost = Kokkos::create_mirror_view(columnDensityRight); + Kokkos::deep_copy(columnDensityLeftHost,columnDensityLeft); + Kokkos::deep_copy(columnDensityRightHost,columnDensityRight); + + errMax = 0.0; + + for(int k = data.beg[KDIR]; k < data.end[KDIR] ; k++) { + for(int j = data.beg[JDIR]; j < data.end[JDIR] ; j++) { + for(int i = data.beg[IDIR]; i < data.end[IDIR] ; i++) { + real err = std::fabs(columnDensityLeftHost(k,j,i)-d.xr[KDIR](k)); + if(err>errMax) errMax=err; + err = std::fabs(columnDensityRightHost(k,j,i)-(1-d.xl[KDIR](k))); + if(err>errMax) errMax=err; + } + } + } + idfx::cout << "Error on column density in KDIR=" << std::scientific << errMax << std::endl; + if(errMax>1e-14) { + IDEFIX_ERROR("Error above tolerance"); + } + +} + +void InternalBoundary(Fluid * hydro, const real t) { + IdefixArray4D Vc = hydro->Vc; + idefix_for("InternalBoundary",0,hydro->data->np_tot[KDIR], + 0,hydro->data->np_tot[JDIR], + 0,hydro->data->np_tot[IDIR], + KOKKOS_LAMBDA (int k, int j, int i) { + // Cancel any motion that could be happening + Vc(VX1,k,j,i) = 0.0; + Vc(VX2,k,j,i) = 0.0; + Vc(VX3,k,j,i) = 0.0; + }); +} + + +// Initialisation routine. Can be used to allocate +// Arrays or variables which are used later on +Setup::Setup(Input &input, Grid &grid, DataBlock &data, Output &output) { + output.EnrollAnalysis(&Analysis); + data.hydro->EnrollInternalBoundary(&InternalBoundary); + + columnX1Left = new Column(IDIR, 1, &data); + columnX1Right = new Column(IDIR, -1, &data); + columnX2Left = new Column(JDIR, 1, &data); + columnX2Right = new Column(JDIR, -1, &data); + columnX3Left = new Column(KDIR, 1, &data); + columnX3Right = new Column(KDIR, -1, &data); + // Initialise the output file +} + +// This routine initialize the flow +// Note that data is on the device. +// One can therefore define locally +// a datahost and sync it, if needed +void Setup::InitFlow(DataBlock &data) { + // Create a host copy + DataBlockHost d(data); + + + for(int k = 0; k < d.np_tot[KDIR] ; k++) { + for(int j = 0; j < d.np_tot[JDIR] ; j++) { + for(int i = 0; i < d.np_tot[IDIR] ; i++) { + + d.Vc(RHO,k,j,i) = 1.0; + d.Vc(VX1,k,j,i) = ZERO_F; + d.Vc(PRS,k,j,i) = 1.0; + + } + } + } + + // Send it all, if needed + d.SyncToDevice(); +} diff --git a/test/utils/columnDensity/testme.py b/test/utils/columnDensity/testme.py new file mode 100755 index 00000000..6f017d7d --- /dev/null +++ b/test/utils/columnDensity/testme.py @@ -0,0 +1,23 @@ +#!/usr/bin/env python3 + +""" + +@author: glesur +""" +import os +import sys +sys.path.append(os.getenv("IDEFIX_DIR")) +import pytools.idfx_test as tst + +test=tst.idfxTest() + +test.configure() +test.compile() +# this test succeeds if it runs successfully +test.run() + +test.mpi = True +test.configure() +test.compile() +# this test succeeds if it runs successfully +test.run()