From e2f7a9a16bf97e2f76245757bbf91cdba188a120 Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Sun, 12 Jan 2025 18:07:36 +0100 Subject: [PATCH 1/8] add a column Density example --- test/utils/columnDensity/column.cpp | 215 +++++++++++++++++++++++ test/utils/columnDensity/column.hpp | 48 +++++ test/utils/columnDensity/definitions.hpp | 4 + test/utils/columnDensity/idefix.ini | 24 +++ test/utils/columnDensity/setup.cpp | 154 ++++++++++++++++ 5 files changed, 445 insertions(+) create mode 100644 test/utils/columnDensity/column.cpp create mode 100644 test/utils/columnDensity/column.hpp create mode 100644 test/utils/columnDensity/definitions.hpp create mode 100644 test/utils/columnDensity/idefix.ini create mode 100644 test/utils/columnDensity/setup.cpp diff --git a/test/utils/columnDensity/column.cpp b/test/utils/columnDensity/column.cpp new file mode 100644 index 000000000..9ad4c3573 --- /dev/null +++ b/test/utils/columnDensity/column.cpp @@ -0,0 +1,215 @@ +#include "column.hpp" +#include "loop.hpp" + +Column::Column(int dir, int sign, int variable, DataBlock *data) : direction(dir), sign(sign), variable(variable) { + 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) { + 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; + + const int var = this->variable; + 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 + this->arr4D = IdefixArray4D (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(this->arr4D); + break; + case 1: + this->mpi.ExchangeX2(this->arr4D); + break; + case 2: + this->mpi.ExchangeX3(this->arr4D); + break; + } + } + } + #endif + idfx::popRegion(); +} diff --git a/test/utils/columnDensity/column.hpp b/test/utils/columnDensity/column.hpp new file mode 100644 index 000000000..b4cc5b539 --- /dev/null +++ b/test/utils/columnDensity/column.hpp @@ -0,0 +1,48 @@ +#ifndef COLUMN_HPP_ +#define COLUMN_HPP_ + +#include "idefix.hpp" +#include "input.hpp" +#include "output.hpp" +#include "grid.hpp" +#include "dataBlock.hpp" +#include "dataBlockHost.hpp" +#include +#include + + + +class Column { + public: + // Constructor from Setup arguments + Column(int dir, int sign, int variable, DataBlock *); + void ComputeColumn(IdefixArray4D in); + 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 + int variable; // The variable we use to compute the column + 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; + + IdefixArray4D arr4D; // Dummy array 4D used for MPI xchanges + std::array nproc; // 3D size of the MPI cartesian geometry + + #endif +}; + +#endif // COLUMN_HPP_ diff --git a/test/utils/columnDensity/definitions.hpp b/test/utils/columnDensity/definitions.hpp new file mode 100644 index 000000000..3b581bd62 --- /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 000000000..2d994315c --- /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 64 u 1.0 +X3-grid 1 0.0 64 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 000000000..73d98cbae --- /dev/null +++ b/test/utils/columnDensity/setup.cpp @@ -0,0 +1,154 @@ +#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); + + columnX1Left->ComputeColumn(data.hydro->Vc); + columnX1Right->ComputeColumn(data.hydro->Vc); + columnX2Left->ComputeColumn(data.hydro->Vc); + columnX2Right->ComputeColumn(data.hydro->Vc); + columnX3Left->ComputeColumn(data.hydro->Vc); + columnX3Right->ComputeColumn(data.hydro->Vc); + + 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; + + 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, RHO, &data); + columnX1Right = new Column(IDIR, -1, RHO, &data); + columnX2Left = new Column(JDIR, 1, RHO, &data); + columnX2Right = new Column(JDIR, -1, RHO, &data); + columnX3Left = new Column(KDIR, 1, RHO, &data); + columnX3Right = new Column(KDIR, -1, RHO, &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(); +} From 26ffc60820205069aa7c810cbd2375e42a2ac855 Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Sun, 12 Jan 2025 18:15:31 +0100 Subject: [PATCH 2/8] add some documentation --- test/utils/columnDensity/column.hpp | 8 ++++++-- test/utils/columnDensity/setup.cpp | 5 +++++ 2 files changed, 11 insertions(+), 2 deletions(-) diff --git a/test/utils/columnDensity/column.hpp b/test/utils/columnDensity/column.hpp index b4cc5b539..9defe7c6b 100644 --- a/test/utils/columnDensity/column.hpp +++ b/test/utils/columnDensity/column.hpp @@ -11,15 +11,19 @@ #include - +// A class to implement a parralel cumulative sum class Column { public: // Constructor from Setup arguments Column(int dir, int sign, int variable, DataBlock *); + // dir : direction along which the integration is performed + // sign: +1 for an integration from left to right, -1 for an integration from right to left (backwards) + // variable: index of the variable along which we do the integral (since the intput array is 4D) void ComputeColumn(IdefixArray4D in); + // Effectively compute integral from the input array in parameter IdefixArray3D GetColumn() { return (this->ColumnArray); - } + }// Return the column density private: IdefixArray3D ColumnArray; int direction; // The direction along which the column is computed diff --git a/test/utils/columnDensity/setup.cpp b/test/utils/columnDensity/setup.cpp index 73d98cbae..ebc131de7 100644 --- a/test/utils/columnDensity/setup.cpp +++ b/test/utils/columnDensity/setup.cpp @@ -34,6 +34,11 @@ void Analysis(DataBlock & data) { 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++) { From a1e6e674743123cb73eee86d4833faaebc73e743 Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Sun, 12 Jan 2025 18:18:40 +0100 Subject: [PATCH 3/8] add it to the tests --- .github/workflows/idefix-ci-jobs.yml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.github/workflows/idefix-ci-jobs.yml b/.github/workflows/idefix-ci-jobs.yml index c0e622966..dfb89e0fb 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 From 9c58596070646975bb7c3c3295904aa80ee34c1f Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Sun, 12 Jan 2025 18:36:28 +0100 Subject: [PATCH 4/8] imrove documentation --- src/utils/CMakeLists.txt | 2 + .../columnDensity => src/utils}/column.cpp | 30 ++++++-- src/utils/column.hpp | 69 +++++++++++++++++++ test/utils/columnDensity/column.hpp | 52 -------------- 4 files changed, 96 insertions(+), 57 deletions(-) rename {test/utils/columnDensity => src/utils}/column.cpp (89%) create mode 100644 src/utils/column.hpp delete mode 100644 test/utils/columnDensity/column.hpp diff --git a/src/utils/CMakeLists.txt b/src/utils/CMakeLists.txt index b36f8a08e..afb4f57d4 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/test/utils/columnDensity/column.cpp b/src/utils/column.cpp similarity index 89% rename from test/utils/columnDensity/column.cpp rename to src/utils/column.cpp index 9ad4c3573..838f1d9ce 100644 --- a/test/utils/columnDensity/column.cpp +++ b/src/utils/column.cpp @@ -1,7 +1,24 @@ +// *********************************************************************************** +// 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 "loop.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, int variable, DataBlock *data) : direction(dir), sign(sign), variable(variable) { +Column::Column(int dir, int sign, int variable, DataBlock *data) + : direction(dir), sign(sign), variable(variable) { idfx::pushRegion("Column::Column"); this->np_tot = data->np_tot; this->np_int = data->np_int; @@ -157,19 +174,22 @@ void Column::ComputeColumn(IdefixArray4D in) { 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))); + 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))); + 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))); + 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))); }); } diff --git a/src/utils/column.hpp b/src/utils/column.hpp new file mode 100644 index 000000000..feb0a8537 --- /dev/null +++ b/src/utils/column.hpp @@ -0,0 +1,69 @@ +// *********************************************************************************** +// 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) + /// @param variable: index of the variable along which we do the integral (since the + /// intput array is 4D) + /////////////////////////////////////////////////////////////////////////////////// + Column(int dir, int sign, int variable, DataBlock *); + + /////////////////////////////////////////////////////////////////////////////////// + /// @brief Effectively compute integral from the input array in argument + /////////////////////////////////////////////////////////////////////////////////// + void ComputeColumn(IdefixArray4D 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 + int variable; // The variable we use to compute the column + 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; + + IdefixArray4D arr4D; // Dummy array 4D used for MPI xchanges + std::array nproc; // 3D size of the MPI cartesian geometry + + #endif +}; + +#endif // UTILS_COLUMN_HPP_ diff --git a/test/utils/columnDensity/column.hpp b/test/utils/columnDensity/column.hpp deleted file mode 100644 index 9defe7c6b..000000000 --- a/test/utils/columnDensity/column.hpp +++ /dev/null @@ -1,52 +0,0 @@ -#ifndef COLUMN_HPP_ -#define COLUMN_HPP_ - -#include "idefix.hpp" -#include "input.hpp" -#include "output.hpp" -#include "grid.hpp" -#include "dataBlock.hpp" -#include "dataBlockHost.hpp" -#include -#include - - -// A class to implement a parralel cumulative sum -class Column { - public: - // Constructor from Setup arguments - Column(int dir, int sign, int variable, DataBlock *); - // dir : direction along which the integration is performed - // sign: +1 for an integration from left to right, -1 for an integration from right to left (backwards) - // variable: index of the variable along which we do the integral (since the intput array is 4D) - void ComputeColumn(IdefixArray4D in); - // Effectively compute integral from the input array in parameter - IdefixArray3D GetColumn() { - return (this->ColumnArray); - }// Return the column density - 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 - int variable; // The variable we use to compute the column - 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; - - IdefixArray4D arr4D; // Dummy array 4D used for MPI xchanges - std::array nproc; // 3D size of the MPI cartesian geometry - - #endif -}; - -#endif // COLUMN_HPP_ From 6283fd5e2997d3502f81b43cfc9110c76b3216f7 Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Sun, 12 Jan 2025 18:50:44 +0100 Subject: [PATCH 5/8] add doc link for column density --- doc/source/faq.rst | 3 +++ 1 file changed, 3 insertions(+) diff --git a/doc/source/faq.rst b/doc/source/faq.rst index 75a07a68b..4c7bc056b 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`) From 8277ed0fdff7114f2c77ff8b30abfa30972c9cab Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Sun, 12 Jan 2025 20:16:30 +0100 Subject: [PATCH 6/8] missing testme file --- test/utils/columnDensity/testme.py | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) create mode 100755 test/utils/columnDensity/testme.py diff --git a/test/utils/columnDensity/testme.py b/test/utils/columnDensity/testme.py new file mode 100755 index 000000000..6f017d7d8 --- /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() From 31ae1b8a47168fcf09cf1f4be6d6794d5a19b33e Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Tue, 14 Jan 2025 15:13:38 +0100 Subject: [PATCH 7/8] add a 3D interface for column density computation --- src/utils/column.cpp | 21 +++++++++++------- src/utils/column.hpp | 15 ++++++++----- test/utils/columnDensity/setup.cpp | 34 +++++++++++++++++++----------- 3 files changed, 45 insertions(+), 25 deletions(-) diff --git a/src/utils/column.cpp b/src/utils/column.cpp index 838f1d9ce..a774d4927 100644 --- a/src/utils/column.cpp +++ b/src/utils/column.cpp @@ -17,8 +17,8 @@ #include "dataBlock.hpp" #include "dataBlockHost.hpp" -Column::Column(int dir, int sign, int variable, DataBlock *data) - : direction(dir), sign(sign), variable(variable) { +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; @@ -62,7 +62,7 @@ Column::Column(int dir, int sign, int variable, DataBlock *data) idfx::popRegion(); } -void Column::ComputeColumn(IdefixArray4D in) { +void Column::ComputeColumn(IdefixArray4D in, const int var) { idfx::pushRegion("Column::ComputeColumn"); const int nk = np_int[KDIR]; const int nj = np_int[JDIR]; @@ -80,7 +80,6 @@ void Column::ComputeColumn(IdefixArray4D in) { auto A = this->Area; auto localSum = this->localSum; - const int var = this->variable; if(direction==IDIR) { // Inspired from loop.hpp Kokkos::parallel_for("ColumnX1", team_policy (nk*nj, Kokkos::AUTO), @@ -210,7 +209,7 @@ void Column::ComputeColumn(IdefixArray4D in) { // density in the ghost zones are coherent #ifdef WITH_MPI // Create a 4D array that contains our column data - this->arr4D = IdefixArray4D (column.data(), 1, this->np_tot[KDIR], + IdefixArray4D arr4D(column.data(), 1, this->np_tot[KDIR], this->np_tot[JDIR], this->np_tot[IDIR]); @@ -219,13 +218,13 @@ void Column::ComputeColumn(IdefixArray4D in) { if(nproc[dir]>1) { switch(dir) { case 0: - this->mpi.ExchangeX1(this->arr4D); + this->mpi.ExchangeX1(arr4D); break; case 1: - this->mpi.ExchangeX2(this->arr4D); + this->mpi.ExchangeX2(arr4D); break; case 2: - this->mpi.ExchangeX3(this->arr4D); + this->mpi.ExchangeX3(arr4D); break; } } @@ -233,3 +232,9 @@ void Column::ComputeColumn(IdefixArray4D in) { #endif idfx::popRegion(); } + +void Column::ComputeColumn(IdefixArray3D in) { + // 4D alias + IdefixArray4D arr4D(in.data(), 1, in.extent(2), in.extent(1), in.extent(0)); + return this->ComputeColumn(arr4D,0); +} diff --git a/src/utils/column.hpp b/src/utils/column.hpp index feb0a8537..43383cbe6 100644 --- a/src/utils/column.hpp +++ b/src/utils/column.hpp @@ -24,15 +24,22 @@ class Column { /// @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) /////////////////////////////////////////////////////////////////////////////////// - Column(int dir, int sign, int variable, DataBlock *); + void ComputeColumn(IdefixArray4D in, int variable); - /////////////////////////////////////////////////////////////////////////////////// + /////////////////////////////////////////////////////////////////////////////////// /// @brief Effectively compute integral from the input array in argument + /// @param in: 3D input array /////////////////////////////////////////////////////////////////////////////////// - void ComputeColumn(IdefixArray4D in); + void ComputeColumn(IdefixArray3D in); /////////////////////////////////////////////////////////////////////////////////// /// @brief Get a reference to the computed column density array @@ -45,7 +52,6 @@ class Column { IdefixArray3D ColumnArray; int direction; // The direction along which the column is computed int sign; // whether we integrate from the left or from the right - int variable; // The variable we use to compute the column std::array np_tot; std::array np_int; std::array beg; @@ -60,7 +66,6 @@ class Column { int MPIrank; int MPIsize; - IdefixArray4D arr4D; // Dummy array 4D used for MPI xchanges std::array nproc; // 3D size of the MPI cartesian geometry #endif diff --git a/test/utils/columnDensity/setup.cpp b/test/utils/columnDensity/setup.cpp index ebc131de7..b27b6b1f3 100644 --- a/test/utils/columnDensity/setup.cpp +++ b/test/utils/columnDensity/setup.cpp @@ -15,12 +15,22 @@ void Analysis(DataBlock & data) { DataBlockHost d(data); - columnX1Left->ComputeColumn(data.hydro->Vc); - columnX1Right->ComputeColumn(data.hydro->Vc); - columnX2Left->ComputeColumn(data.hydro->Vc); - columnX2Right->ComputeColumn(data.hydro->Vc); - columnX3Left->ComputeColumn(data.hydro->Vc); - columnX3Right->ComputeColumn(data.hydro->Vc); + // 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",0,data.np_tot[KDIR],0,data.np_tot[JDIR],0,data.np_tot[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; @@ -124,12 +134,12 @@ Setup::Setup(Input &input, Grid &grid, DataBlock &data, Output &output) { output.EnrollAnalysis(&Analysis); data.hydro->EnrollInternalBoundary(&InternalBoundary); - columnX1Left = new Column(IDIR, 1, RHO, &data); - columnX1Right = new Column(IDIR, -1, RHO, &data); - columnX2Left = new Column(JDIR, 1, RHO, &data); - columnX2Right = new Column(JDIR, -1, RHO, &data); - columnX3Left = new Column(KDIR, 1, RHO, &data); - columnX3Right = new Column(KDIR, -1, RHO, &data); + 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 } From 7c10ec052d60128aea8d4dd90ee32c8a13f9b06b Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Wed, 15 Jan 2025 18:06:53 +0100 Subject: [PATCH 8/8] fix dimensions for alias --- src/utils/column.cpp | 2 +- test/utils/columnDensity/idefix.ini | 4 ++-- test/utils/columnDensity/setup.cpp | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/utils/column.cpp b/src/utils/column.cpp index a774d4927..34d55f56b 100644 --- a/src/utils/column.cpp +++ b/src/utils/column.cpp @@ -235,6 +235,6 @@ void Column::ComputeColumn(IdefixArray4D in, const int var) { void Column::ComputeColumn(IdefixArray3D in) { // 4D alias - IdefixArray4D arr4D(in.data(), 1, in.extent(2), in.extent(1), in.extent(0)); + IdefixArray4D arr4D(in.data(), 1, in.extent(0), in.extent(1), in.extent(2)); return this->ComputeColumn(arr4D,0); } diff --git a/test/utils/columnDensity/idefix.ini b/test/utils/columnDensity/idefix.ini index 2d994315c..5e80c8bce 100644 --- a/test/utils/columnDensity/idefix.ini +++ b/test/utils/columnDensity/idefix.ini @@ -1,7 +1,7 @@ [Grid] X1-grid 1 0.0 64 u 1.0 -X2-grid 1 0.0 64 u 1.0 -X3-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 diff --git a/test/utils/columnDensity/setup.cpp b/test/utils/columnDensity/setup.cpp index b27b6b1f3..fe7c57e69 100644 --- a/test/utils/columnDensity/setup.cpp +++ b/test/utils/columnDensity/setup.cpp @@ -23,7 +23,7 @@ void Analysis(DataBlock & data) { IdefixArray3D rho("rho",data.np_tot[KDIR],data.np_tot[JDIR],data.np_tot[IDIR]); auto Vc = data.hydro->Vc; - idefix_for("init rho",0,data.np_tot[KDIR],0,data.np_tot[JDIR],0,data.np_tot[IDIR], + 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); });