diff --git a/test/MHD/OrszagTang/idefix-hll.ini b/test/MHD/OrszagTang/idefix-hll.ini index 4f13fd972..707f35605 100644 --- a/test/MHD/OrszagTang/idefix-hll.ini +++ b/test/MHD/OrszagTang/idefix-hll.ini @@ -18,6 +18,7 @@ X2-beg periodic X2-end periodic [Output] -vtk 0.5 -dmp 0.5 -log 100 +analysis 0.1 +vtk 0.5 +dmp 0.5 +log 100 diff --git a/test/MHD/OrszagTang/idefix-hlld-arithmetic.ini b/test/MHD/OrszagTang/idefix-hlld-arithmetic.ini index 8a656d7b3..0c62cfcc5 100644 --- a/test/MHD/OrszagTang/idefix-hlld-arithmetic.ini +++ b/test/MHD/OrszagTang/idefix-hlld-arithmetic.ini @@ -19,6 +19,7 @@ X2-beg periodic X2-end periodic [Output] -vtk 0.5 -dmp 0.5 -log 100 +analysis 0.1 +vtk 0.5 +dmp 0.5 +log 100 diff --git a/test/MHD/OrszagTang/idefix-hlld-hll.ini b/test/MHD/OrszagTang/idefix-hlld-hll.ini index b396e9917..26ae64973 100644 --- a/test/MHD/OrszagTang/idefix-hlld-hll.ini +++ b/test/MHD/OrszagTang/idefix-hlld-hll.ini @@ -19,6 +19,7 @@ X2-beg periodic X2-end periodic [Output] -vtk 0.5 -dmp 0.5 -log 100 +analysis 0.1 +vtk 0.5 +dmp 0.5 +log 100 diff --git a/test/MHD/OrszagTang/idefix-hlld-hlld.ini b/test/MHD/OrszagTang/idefix-hlld-hlld.ini index 5e158b939..8f82d6f8e 100644 --- a/test/MHD/OrszagTang/idefix-hlld-hlld.ini +++ b/test/MHD/OrszagTang/idefix-hlld-hlld.ini @@ -19,6 +19,7 @@ X2-beg periodic X2-end periodic [Output] -vtk 0.5 -dmp 0.5 -log 100 +analysis 0.1 +vtk 0.5 +dmp 0.5 +log 100 diff --git a/test/MHD/OrszagTang/idefix-hlld-uct0.ini b/test/MHD/OrszagTang/idefix-hlld-uct0.ini index 23646b96f..b0d637d74 100644 --- a/test/MHD/OrszagTang/idefix-hlld-uct0.ini +++ b/test/MHD/OrszagTang/idefix-hlld-uct0.ini @@ -19,6 +19,7 @@ X2-beg periodic X2-end periodic [Output] -vtk 0.5 -dmp 0.5 -log 100 +analysis 0.1 +vtk 0.5 +dmp 0.5 +log 100 diff --git a/test/MHD/OrszagTang/idefix-hlld.ini b/test/MHD/OrszagTang/idefix-hlld.ini index cf6d5cbfd..449d941f4 100644 --- a/test/MHD/OrszagTang/idefix-hlld.ini +++ b/test/MHD/OrszagTang/idefix-hlld.ini @@ -18,6 +18,7 @@ X2-beg periodic X2-end periodic [Output] -vtk 0.5 -dmp 0.5 -log 100 +analysis 0.1 +vtk 0.5 +dmp 0.5 +log 100 diff --git a/test/MHD/OrszagTang/idefix-tvdlf.ini b/test/MHD/OrszagTang/idefix-tvdlf.ini index ec9ffe85e..62cf6d280 100644 --- a/test/MHD/OrszagTang/idefix-tvdlf.ini +++ b/test/MHD/OrszagTang/idefix-tvdlf.ini @@ -18,6 +18,7 @@ X2-beg periodic X2-end periodic [Output] -vtk 0.5 -dmp 0.5 -log 100 +analysis 0.1 +vtk 0.5 +dmp 0.5 +log 100 diff --git a/test/MHD/OrszagTang/idefix.ini b/test/MHD/OrszagTang/idefix.ini index 0e7421843..2ecc050ec 100644 --- a/test/MHD/OrszagTang/idefix.ini +++ b/test/MHD/OrszagTang/idefix.ini @@ -18,6 +18,7 @@ X2-beg periodic X2-end periodic [Output] -vtk 0.5 -dmp 0.5 -log 100 +analysis 0.1 +vtk 0.5 +dmp 0.5 +log 100 diff --git a/test/MHD/OrszagTang/setup.cpp b/test/MHD/OrszagTang/setup.cpp index 50c981c18..b337ff6bf 100644 --- a/test/MHD/OrszagTang/setup.cpp +++ b/test/MHD/OrszagTang/setup.cpp @@ -9,13 +9,53 @@ generators on different architectures. **/ /*********************************************/ - +void CheckConservation(DataBlock &data) { + static real firstCall{true}; + static std::array consArray; + data.hydro->ConvertPrimToCons(); + auto Uc = data.hydro->Uc; + auto dV = data.dV; + //idfx::cout << "Analysis: checking conserved quantities..." << std::endl; + #ifdef SINGLE_PRECISION + const real threshold = 1e-4; + #else + const real threshold = 1e-13; + #endif + for(int nv = 0 ; nv < DefaultPhysics::nvar ; nv++) { + real cons = 0; + idefix_reduce("Conserved quantity reduction", + 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, real &c) { + c += dV(k,j,i)*Uc(nv,k,j,i); + }, + Kokkos::Sum(cons)); + #ifdef WITH_MPI + MPI_Allreduce(MPI_IN_PLACE, &cons, 1, realMPI, MPI_SUM, MPI_COMM_WORLD); + #endif + if(firstCall) { + consArray[nv] = cons; + } else { + real err = std::fabs((consArray[nv]-cons)); + if(err>threshold) { + std::stringstream str; + str << "Quantity " << data.hydro->VcName[nv] << " is not conserved" << std::endl; + std::cout << "Error on " << data.hydro->VcName[nv] << " is " << err << std::endl; + std::cout << "Original=" << consArray[nv] << " New=" << cons << std::endl; + IDEFIX_ERROR(str); + } + } + } + firstCall=false; + //idfx::cout << "Analysis: done." << std::endl; +} // Default constructor // 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(&CheckConservation); } // This routine initialize the flow diff --git a/test/MHD/OrszagTang3D/idefix.ini b/test/MHD/OrszagTang3D/idefix.ini index 1021bd335..5551b23f8 100644 --- a/test/MHD/OrszagTang3D/idefix.ini +++ b/test/MHD/OrszagTang3D/idefix.ini @@ -22,6 +22,7 @@ X3-beg periodic X3-end periodic [Output] -vtk 0.2 -dmp 0.2 -log 10 +vtk 0.2 +dmp 0.2 +log 10 +analysis 0.1 diff --git a/test/MHD/OrszagTang3D/setup.cpp b/test/MHD/OrszagTang3D/setup.cpp index a6927a132..abbc78534 100644 --- a/test/MHD/OrszagTang3D/setup.cpp +++ b/test/MHD/OrszagTang3D/setup.cpp @@ -11,11 +11,53 @@ generators on different architectures. Output* myOutput; int outnum; -// Analysis function + +void CheckConservation(DataBlock &data) { + static real firstCall{true}; + static std::array consArray; + data.hydro->ConvertPrimToCons(); + auto Uc = data.hydro->Uc; + auto dV = data.dV; + //idfx::cout << "Analysis: checking conserved quantities..." << std::endl; + #ifdef SINGLE_PRECISION + const real threshold = 2e-4; + #else + const real threshold = 1e-13; + #endif + for(int nv = 0 ; nv < DefaultPhysics::nvar ; nv++) { + real cons = 0; + idefix_reduce("Conserved quantity reduction", + 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, real &c) { + c += dV(k,j,i)*Uc(nv,k,j,i); + }, + Kokkos::Sum(cons)); + #ifdef WITH_MPI + MPI_Allreduce(MPI_IN_PLACE, &cons, 1, realMPI, MPI_SUM, MPI_COMM_WORLD); + #endif + if(firstCall) { + consArray[nv] = cons; + } else { + real err = std::fabs((consArray[nv]-cons)); + if(err>threshold) { + std::stringstream str; + str << "Quantity " << data.hydro->VcName[nv] << " is not conserved" << std::endl; + std::cout << "Error on " << data.hydro->VcName[nv] << " is " << err << std::endl; + std::cout << "Original=" << consArray[nv] << " New=" << cons << std::endl; + IDEFIX_ERROR(str); + } + } + } + firstCall=false; + //idfx::cout << "Analysis: done." << std::endl; +} // 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(&CheckConservation); } // This routine initialize the flow