From 58bc3f1f892161469e1761e6fae56aec717c6f8b Mon Sep 17 00:00:00 2001 From: u43070nw Date: Fri, 1 Mar 2024 15:16:58 +0000 Subject: [PATCH 1/3] Change to HeFFTe & FFTW3 backend --- CMakeLists.txt | 14 +- app/CMakeLists.txt | 4 + app/OT/BEqn.H | 12 +- app/OT_Poisson/OT.cpp | 19 ++- app/TG/TG.cpp | 2 +- app/TG_Piosson/RKLoop.H | 1 + app/TG_Piosson/TG.cpp | 8 +- app/shearLayer/TEqn.H | 38 +++++ app/shearLayer/UEqn.H | 51 +++++++ app/shearLayer/hydrostaticP.H | 39 +++++ app/shearLayer/layerAverage.H | 52 +++++++ app/shearLayer/pEqn.H | 4 + app/shearLayer/rhoEqn.H | 38 +++++ app/shearLayer/shearLayer.cpp | 195 +++++++++++++++++++++++++ app/shearLayerWC/TEqn.H | 38 +++++ app/shearLayerWC/UEqn.H | 42 ++++++ app/shearLayerWC/hydrostaticP.H | 39 +++++ app/shearLayerWC/layerAverage.H | 52 +++++++ app/shearLayerWC/pEqn.H | 4 + app/shearLayerWC/rhoEqn.H | 38 +++++ app/shearLayerWC/shearLayer.cpp | 247 ++++++++++++++++++++++++++++++++ examples/OT/UBC | 36 ++--- examples/shearLayer/TBC | 6 + examples/shearLayer/UBC | 18 +++ examples/shearLayer/mesh | 4 + examples/shearLayer/pBC | 6 + examples/shearLayer/pHydBC | 6 + findFFTMPI/FindFFTMPI.cmake | 64 --------- findHEFFTE/findHEFFTE.cmake | 27 ++++ src/BC/BC.h | 2 +- src/BC/dirichletBC.hxx | 19 ++- src/BC/neumannBC.hxx | 58 +++++--- src/Mesh/Mesh.cpp | 5 +- src/Mesh/Mesh.h | 4 +- src/parallelCom/parallelCom.cpp | 26 +++- src/poisson/poisson.cpp | 125 +++++++--------- src/poisson/poisson.h | 23 +-- 37 files changed, 1147 insertions(+), 219 deletions(-) create mode 100644 app/shearLayer/TEqn.H create mode 100644 app/shearLayer/UEqn.H create mode 100644 app/shearLayer/hydrostaticP.H create mode 100644 app/shearLayer/layerAverage.H create mode 100644 app/shearLayer/pEqn.H create mode 100644 app/shearLayer/rhoEqn.H create mode 100644 app/shearLayer/shearLayer.cpp create mode 100644 app/shearLayerWC/TEqn.H create mode 100644 app/shearLayerWC/UEqn.H create mode 100644 app/shearLayerWC/hydrostaticP.H create mode 100644 app/shearLayerWC/layerAverage.H create mode 100644 app/shearLayerWC/pEqn.H create mode 100644 app/shearLayerWC/rhoEqn.H create mode 100644 app/shearLayerWC/shearLayer.cpp create mode 100644 examples/shearLayer/TBC create mode 100644 examples/shearLayer/UBC create mode 100644 examples/shearLayer/mesh create mode 100644 examples/shearLayer/pBC create mode 100644 examples/shearLayer/pHydBC delete mode 100644 findFFTMPI/FindFFTMPI.cmake create mode 100644 findHEFFTE/findHEFFTE.cmake diff --git a/CMakeLists.txt b/CMakeLists.txt index 359ebc7..cb11345 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -20,7 +20,7 @@ endif() set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/cmake/Modules/") set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${PROJECT_SOURCE_DIR}/findFFTW") -set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${PROJECT_SOURCE_DIR}/findFFTMPI") +set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${PROJECT_SOURCE_DIR}/findHEFFTE") set(LIBS) set(CMAKE_CXX_STANDARD 20) # C++20... @@ -38,7 +38,7 @@ include_directories(${Boost_INCLUDE_DIRS}) option(USE_MPI "Use MPI" ON) option(USE_OMP "Use OpenMP" OFF) -option(USE_FFTMPI "Use FFTMPI" ON) +option(USE_HEFFTE "Use Heffte" ON) ###mpi if(USE_MPI) @@ -64,13 +64,13 @@ if(USE_OMP) message(SEND_ERROR "OMP not yet supported / experimental.") endif() -if(USE_FFTMPI) +if(USE_HEFFTE) + add_definitions(-DHAVE_HEFFTE) + find_package(Heffte REQUIRED) find_package(FFTW REQUIRED) - find_package(FFTMPI REQUIRED) - add_definitions(-DHAVE_FFTMPI) + include_directories(${HEFFTE_INCLUDES}) include_directories(${FFTW_INCLUDE_DIRS}) - include_directories(${FFTMPI_INCLUDES}) - set(LIBS ${LIBS} "${FFTMPI_LIBRARIES}" "${FFTW_LIBRARIES}" ) + set(LIBS ${LIBS} "${HEFFTE_LIBRARIES}" "${FFTW_LIBRARIES}") endif() option(PROFILE "Use gperftools" OFF) diff --git a/app/CMakeLists.txt b/app/CMakeLists.txt index 3fa2ced..43d9daa 100644 --- a/app/CMakeLists.txt +++ b/app/CMakeLists.txt @@ -13,3 +13,7 @@ install(TARGETS OT DESTINATION ${CMAKE_INSTALL_PREFIX}/bin) add_executable (OT_Poisson ./OT_Poisson/OT.cpp) target_link_libraries (OT_Poisson NOVA ${Boost_LIBRARIES} ${MPI_LIBRARIES} ${FFTMPI_LIBRARIES} ${LIBS}) install(TARGETS OT DESTINATION ${CMAKE_INSTALL_PREFIX}/bin) + +add_executable (shearLayer ./shearLayer/shearLayer.cpp) +target_link_libraries (shearLayer NOVA ${Boost_LIBRARIES} ${MPI_LIBRARIES} ${FFTMPI_LIBRARIES} ${LIBS}) +install(TARGETS shearLayer DESTINATION ${CMAKE_INSTALL_PREFIX}/bin) diff --git a/app/OT/BEqn.H b/app/OT/BEqn.H index 9f6224d..f988aae 100644 --- a/app/OT/BEqn.H +++ b/app/OT/BEqn.H @@ -29,12 +29,12 @@ } - Bstar.correctBoundaryConditions(); - pBEqn.rhs( ex::div(Bstar) ); // /time.dt() - pBEqn.solve(); - pB.correctBoundaryConditions(); - B = Bstar - ex::grad(pB); //*time.dt() - B.correctBoundaryConditions(); +// Bstar.correctBoundaryConditions(); +// pBEqn.rhs( ex::div(Bstar) ); // /time.dt() +// pBEqn.solve(); +// pB.correctBoundaryConditions(); +// B = Bstar - ex::grad(pB); //*time.dt() +// B.correctBoundaryConditions(); if( time.writelog() ) { diff --git a/app/OT_Poisson/OT.cpp b/app/OT_Poisson/OT.cpp index 1fe3e67..9943f45 100644 --- a/app/OT_Poisson/OT.cpp +++ b/app/OT_Poisson/OT.cpp @@ -23,7 +23,7 @@ int main(int argc, char* argv[]) { settings::process( argc, argv ); - Time time( 0.001, 8, 50 ); //args: dt, endT, write interval / steps + Time time( 0.005, 8.01, 20 ); //args: dt, endT, write interval / steps const scalar pi = tools::pi; parallelCom::decompose( settings::zoneName()+"/"+"mesh" ); @@ -51,6 +51,9 @@ int main(int argc, char* argv[]) poisson pEqn(p_ptr); poisson pBEqn(pB_ptr); + std::default_random_engine eng( parallelCom::myProcNo() ); + std::uniform_real_distribution dist(-1.0,1.0); + //initial conditions for( int i=settings::m()/2; i dTt( mesh ); + std::shared_ptr > dT( dTt() ); + + reuseTmp Tct( mesh ); + std::shared_ptr > Tc( Tct() ); + + reuseTmp TOldt( mesh ); + std::shared_ptr > TOld( TOldt() ); + + *Tc = T; + *TOld = T; + + for( int rk=0; rk<4; rk++ ) + { + *dT = + time.dt() * + ( + -ex::div( U, T ) + +ex::laplacian( ka, T ) + ); + +// drhoT->correctBoundaryConditions(); + + *Tc = Tc+RK4[rk]*dT; + + if( rk<3 ) + { + T = TOld + (int(rk/2)+1)/2.0 * (dT); + } + else + { + T=Tc; + } + + T.correctBoundaryConditions(); + } +} diff --git a/app/shearLayer/UEqn.H b/app/shearLayer/UEqn.H new file mode 100644 index 0000000..9c77743 --- /dev/null +++ b/app/shearLayer/UEqn.H @@ -0,0 +1,51 @@ +{ +// std::array >, 4> dU; +// for( int i=0; i<4; i++ ) +// { + reuseTmp dUt( mesh ); + std::shared_ptr > dU( dUt() ); +// dU[i]= dUt(); +// } + + reuseTmp UOldt( mesh ); + std::shared_ptr > UOld( UOldt() ); + + *UOld = U; + +// std::shared_ptr > gradU( fdc::grad( U ) ); +// gradU->correctBoundaryConditions(); + + for( int rk=0; rk<4; rk++ ) + { + *dU = + time.dt() * + ( + -ex::div( U, U ) + +ex::laplacian( mu, U ) +// +mu*ex::div( ( dev2(transpose(gradU) ) ) ) + -( (rho0*T)*g ) + ); + + Ustar = UOld; +// for( int j=0; j<=rk; j++ ) +// { +// Ustar = Ustar+settings::RK4[rk]*dU; + Ustar = Ustar+RK4[rk]*dU; +// } + + Ustar.correctBoundaryConditions(); +// pEqn.rhs(ex::div(Ustar) / time.dt() / settings::RK4c()[rk]); +// pEqn.solve(); + p.correctBoundaryConditions(); +// U = Ustar - settings::RK4c()[rk]*time.dt() * ex::grad(p); //(p_rgh?) + + U.correctBoundaryConditions(); + } + + + + + + + +} diff --git a/app/shearLayer/hydrostaticP.H b/app/shearLayer/hydrostaticP.H new file mode 100644 index 0000000..9f636e3 --- /dev/null +++ b/app/shearLayer/hydrostaticP.H @@ -0,0 +1,39 @@ +{ + for( int i=settings::m()/2; i cumTot( parallelCom::worldSize(), scalar(0.0) ); + cumTot[parallelCom::myProcNo()] = pHydc; + + all_reduce( cumTot, std::plus() ); + + pHydc = 0.0; + double aboveProcs=0.0; + + for( int pk=parallelCom::k()+1; pk=settings::m()/2; k-- ) + { + pHydc += T(i, j, k)*mesh.dz()*g[2]*rho0; + pHyd(i, j, k) = pHydc; // + aboveProcs; + } + + } + } +//add field and implement null boundary condition +// pHyd.correctBoundaryConditions(); + + p += pHyd; +} + diff --git a/app/shearLayer/layerAverage.H b/app/shearLayer/layerAverage.H new file mode 100644 index 0000000..0e7087e --- /dev/null +++ b/app/shearLayer/layerAverage.H @@ -0,0 +1,52 @@ +{ + for( int k=0; k proc_avgU( parallelCom::worldSize(), 0.0 ); + std::vector proc_avgT( parallelCom::worldSize(), 0.0 ); + std::vector proc_N( parallelCom::worldSize(), 0 ); + + for( int j=settings::m()/2; j()); + all_reduce(proc_avgT, std::plus()); + all_reduce(proc_N, std::plus()); + + vector avgU(0,0,0); + scalar avgT=0.0; + int N=0; + + for( int proci=0; procicorrectBoundaryConditions(); + //Tp->correctBoundaryConditions(); + } +} + diff --git a/app/shearLayer/pEqn.H b/app/shearLayer/pEqn.H new file mode 100644 index 0000000..63303fe --- /dev/null +++ b/app/shearLayer/pEqn.H @@ -0,0 +1,4 @@ +{ + p = (c*c/2.0)*rho0*rho0; +} + diff --git a/app/shearLayer/rhoEqn.H b/app/shearLayer/rhoEqn.H new file mode 100644 index 0000000..4d4334c --- /dev/null +++ b/app/shearLayer/rhoEqn.H @@ -0,0 +1,38 @@ +{ + reuseTmp drhot( mesh ); + std::shared_ptr > drho( drhot() ); + + reuseTmp rhoct( mesh ); + std::shared_ptr > rhoc( rhoct() ); + + reuseTmp rhoOldt( mesh ); + std::shared_ptr > rhoOld( rhoOldt() ); + + *rhoOld = rho; + *rhoc = rho; + + for( int rk=0; rk<4; rk++ ) + { + *drho = + time.dt() * + ( + -ex::div( rho*U ) + ); + + //drho->correctBoundaryConditions(); + + *rhoc = rhoc+RK4[rk]*drho; + + if( rk<3 ) + { + rho = rhoOld + (int(rk/2)+1)/2.0 * drho; + } + else + { + rho=rhoc; + } + + rho.correctBoundaryConditions(); + } + +} diff --git a/app/shearLayer/shearLayer.cpp b/app/shearLayer/shearLayer.cpp new file mode 100644 index 0000000..6950223 --- /dev/null +++ b/app/shearLayer/shearLayer.cpp @@ -0,0 +1,195 @@ +#include "Field/Field.h" +#include +#include "Types/vector/vector.h" +#include "ex/grad/grad.h" +#include "ex/div/div.h" +#include "ex/curl/curl.h" +#include "ex/laplacian/laplacian.h" +#include "Time/Time.h" +#include "parallelCom/parallelCom.h" +#include "settings/settings.h" +#include "BC/BC.h" +#include +#include "Tools/Tools.h" +#include "poisson/poisson.h" + +#include + +#include +#include +#include + +int main(int argc, char* argv[]) +{ + settings::process( argc, argv ); + Time time( 0.1, 1000000, 100 ); //args: dt, endT, write interval / steps + + parallelCom::decompose( settings::zoneName()+"/"+"mesh" ); + + Mesh mesh( settings::zoneName()+"/"+"mesh", time ); + + mesh.write(settings::zoneName()+"/data"); + + scalar Mach=0.2; + scalar Re0=1967.0; + scalar Ri0=0.08; + scalar Pr=7; + + scalar rho0=999.57; + scalar c=0.0417; + + scalar u0=c*Mach; + scalar h0=0.2;//Re0*mu/rho0/u0; + scalar t0=0.1;//-Ri0*u0*u0/g[2]/h0; + scalar mu = rho0*u0*h0/Re0; + scalar ka = mu/Pr; + vector g( 0, 0, -Ri0*u0*u0/t0/h0 ); + + Field U( mesh, vector(1,0,0), "U" ); + Field T( mesh, 0, "T" ); + Field pHyd( mesh,0, "pHyd" ); + Field Ustar( mesh, vector(0,0,0), "U" ); +//------------------------------------------- +// time.dt() = 0.5*std::min(mesh.dx(), std::min( mesh.dy(), mesh.dz() ) )/(u0/Mach); +// +// if( parallelCom::master() ) +// { +// std::cout<<"Speed of sound: "< dist(-1.0,1.0); + + std::shared_ptr > p_ptr( std::make_shared >( mesh, 0, "p" ) ); + auto& p = (*p_ptr); + + poisson pEqn(p_ptr); + + { + #include "pEqn.H" + //initial conditions + for( int i=settings::m()/2; i upt( mesh ); +// std::shared_ptr > up( upt() ); +// reuseTmp Tpt( mesh ); +// std::shared_ptr > Tp( Tpt() ); +// #include "layerAverage.H" +// +// reuseTmp epst( mesh ); +// std::shared_ptr > eps( epst() ); +// *eps = 2.0*( ex::grad( up ) && ex::grad( up ) )*mu/rho0; +// (*eps).write("e"); +// +// reuseTmp xit( mesh ); +// std::shared_ptr > xi( xit() ); +// *xi = 2.0*( ex::grad( Tp ) & ex::grad( Tp ) )*ka/rho0; +// (*xi).write("x"); +// +// reuseTmp vortt( mesh ); +// std::shared_ptr > vort( vortt() ); +// *vort = ex::curl( U ); +// (*vort).write("v"); +// +// reuseTmp gpt( mesh ); +// std::shared_ptr > gp( gpt() ); +// *gp = ex::grad(p); +// (*gp).write("z"); +// +// mesh.writeSOS(); +// } + } + + + std::cout<< timer.elapsed().wall / 1e9 < drhoTt( mesh ); + std::shared_ptr > drhoT( drhoTt() ); + + reuseTmp rhoTct( mesh ); + std::shared_ptr > rhoTc( rhoTct() ); + + reuseTmp rhoTOldt( mesh ); + std::shared_ptr > rhoTOld( rhoTOldt() ); + + *rhoTc = rho*T; + *rhoTOld = rho*T; + + for( int rk=0; rk<4; rk++ ) + { + *drhoT = + time.dt() * + ( + -ex::div( rho*U, T ) + +ex::laplacian( ka, T ) + ); + +// drhoT->correctBoundaryConditions(); + + *rhoTc = rhoTc+RK4[rk]*drhoT; + + if( rk<3 ) + { + T = rhoTOld/rho + (int(rk/2)+1)/2.0 * (drhoT/rho); + } + else + { + T=rhoTc/rho; + } + + T.correctBoundaryConditions(); + } +} diff --git a/app/shearLayerWC/UEqn.H b/app/shearLayerWC/UEqn.H new file mode 100644 index 0000000..025635c --- /dev/null +++ b/app/shearLayerWC/UEqn.H @@ -0,0 +1,42 @@ +{ + reuseTmp drhoUt( mesh ); + std::shared_ptr > drhoU( drhoUt() ); + + reuseTmp rhoUOldt( mesh ); + std::shared_ptr > rhoUOld( rhoUOldt() ); + + reuseTmp rhoUct( mesh ); + std::shared_ptr > rhoUc( rhoUct() ); + + *rhoUc = rho*U; + *rhoUOld = rho*U; + +// std::shared_ptr > gradU( fdc::grad( U ) ); +// gradU->correctBoundaryConditions(); + + for( int rk=0; rk<4; rk++ ) + { + *drhoU = + time.dt() * + ( + -ex::div( rho*U, U ) + -ex::grad( p ) + +ex::laplacian( mu, U ) +// +mu*ex::div( ( dev2(transpose(gradU) ) ) ) + -( (rho*T)*g ) + ); + + *rhoUc = rhoUc+RK4[rk]*drhoU; + + if( rk<3 ) + { + U = rhoUOld/rho + (int(rk/2)+1)/2.0 * (drhoU/rho); + } + else + { + U=rhoUc/rho; + } + + U.correctBoundaryConditions(); + } +} diff --git a/app/shearLayerWC/hydrostaticP.H b/app/shearLayerWC/hydrostaticP.H new file mode 100644 index 0000000..9f636e3 --- /dev/null +++ b/app/shearLayerWC/hydrostaticP.H @@ -0,0 +1,39 @@ +{ + for( int i=settings::m()/2; i cumTot( parallelCom::worldSize(), scalar(0.0) ); + cumTot[parallelCom::myProcNo()] = pHydc; + + all_reduce( cumTot, std::plus() ); + + pHydc = 0.0; + double aboveProcs=0.0; + + for( int pk=parallelCom::k()+1; pk=settings::m()/2; k-- ) + { + pHydc += T(i, j, k)*mesh.dz()*g[2]*rho0; + pHyd(i, j, k) = pHydc; // + aboveProcs; + } + + } + } +//add field and implement null boundary condition +// pHyd.correctBoundaryConditions(); + + p += pHyd; +} + diff --git a/app/shearLayerWC/layerAverage.H b/app/shearLayerWC/layerAverage.H new file mode 100644 index 0000000..0e7087e --- /dev/null +++ b/app/shearLayerWC/layerAverage.H @@ -0,0 +1,52 @@ +{ + for( int k=0; k proc_avgU( parallelCom::worldSize(), 0.0 ); + std::vector proc_avgT( parallelCom::worldSize(), 0.0 ); + std::vector proc_N( parallelCom::worldSize(), 0 ); + + for( int j=settings::m()/2; j()); + all_reduce(proc_avgT, std::plus()); + all_reduce(proc_N, std::plus()); + + vector avgU(0,0,0); + scalar avgT=0.0; + int N=0; + + for( int proci=0; procicorrectBoundaryConditions(); + //Tp->correctBoundaryConditions(); + } +} + diff --git a/app/shearLayerWC/pEqn.H b/app/shearLayerWC/pEqn.H new file mode 100644 index 0000000..cd6ff37 --- /dev/null +++ b/app/shearLayerWC/pEqn.H @@ -0,0 +1,4 @@ +{ + p = (c*c/2.0/rho0)*rho*rho; +} + diff --git a/app/shearLayerWC/rhoEqn.H b/app/shearLayerWC/rhoEqn.H new file mode 100644 index 0000000..4d4334c --- /dev/null +++ b/app/shearLayerWC/rhoEqn.H @@ -0,0 +1,38 @@ +{ + reuseTmp drhot( mesh ); + std::shared_ptr > drho( drhot() ); + + reuseTmp rhoct( mesh ); + std::shared_ptr > rhoc( rhoct() ); + + reuseTmp rhoOldt( mesh ); + std::shared_ptr > rhoOld( rhoOldt() ); + + *rhoOld = rho; + *rhoc = rho; + + for( int rk=0; rk<4; rk++ ) + { + *drho = + time.dt() * + ( + -ex::div( rho*U ) + ); + + //drho->correctBoundaryConditions(); + + *rhoc = rhoc+RK4[rk]*drho; + + if( rk<3 ) + { + rho = rhoOld + (int(rk/2)+1)/2.0 * drho; + } + else + { + rho=rhoc; + } + + rho.correctBoundaryConditions(); + } + +} diff --git a/app/shearLayerWC/shearLayer.cpp b/app/shearLayerWC/shearLayer.cpp new file mode 100644 index 0000000..3c04145 --- /dev/null +++ b/app/shearLayerWC/shearLayer.cpp @@ -0,0 +1,247 @@ +#include "Field/Field.h" +#include +#include "Types/vector/vector.h" +#include "ex/grad/grad.h" +#include "ex/div/div.h" +#include "ex/curl/curl.h" +#include "ex/laplacian/laplacian.h" +#include "Time/Time.h" +#include "parallelCom/parallelCom.h" +#include "settings/settings.h" +#include "BC/BC.h" +#include +#include "Tools/Tools.h" + +#include + +#include +#include + +int main(int argc, char* argv[]) +{ + settings::process( argc, argv ); + Time time( 0.01, 1000000, 100 ); //args: dt, endT, write interval / steps + + parallelCom::decompose( settings::zoneName()+"/"+"mesh" ); + + Mesh mesh( settings::zoneName()+"/"+"mesh", time ); + + mesh.write(settings::zoneName()+"/data"); + + scalar Mach=0.2; + scalar Re0=1967.0; + scalar Ri0=0.08; + scalar Pr=7; + + scalar rho0=999.57; + scalar c=0.0417; + + scalar u0=c*Mach; + scalar h0=0.2;//Re0*mu/rho0/u0; + scalar t0=0.1;//-Ri0*u0*u0/g[2]/h0; + scalar mu = rho0*u0*h0/Re0; + scalar ka = mu/Pr; + vector g( 0, 0, -Ri0*u0*u0/t0/h0 ); + + Field U( mesh, vector(1,0,0), "U" ); + Field T( mesh, 0, "T" ); + Field p( mesh,0, "p" ); + Field pHyd( mesh,0, "pHyd" ); + Field rho( mesh, rho0, "rho" ); + +//------------------------------------------- +// Time time( 5e-9, 1000000, scalar(6*h0/u0) ); +// +// int ni=840*1.5; +// int nj=102*1.5; +// int nk=900*1.5; +// +// parallelCom::decompose( ni, nj, nk ); +// +// Mesh mesh( ni, nj, nk, 14.0*h0, 1.7*h0, 15*h0, time ); +// +// +// time.dt() = 0.5*std::min(mesh.dx(), std::min( mesh.dy(), mesh.dz() ) )/(u0/Mach); +// +// if( parallelCom::master() ) +// { +// std::cout<<"Speed of sound: "< dist(-1.0,1.0); + + +// if( settings::restart() ) +// { +// U.read("U"); +// rho.read("r"); +// T.read("T"); +// p.read("p"); +// +// U.correctBoundaryConditions(); +// p.correctBoundaryConditions(); +// rho.correctBoundaryConditions(); +// T.correctBoundaryConditions(); +// } +// else + { + #include "pEqn.H" + //initial conditions + for( int i=settings::m()/2; i upt( mesh ); +// std::shared_ptr > up( upt() ); +// reuseTmp Tpt( mesh ); +// std::shared_ptr > Tp( Tpt() ); +// #include "layerAverage.H" +// +// reuseTmp epst( mesh ); +// std::shared_ptr > eps( epst() ); +// *eps = 2.0*( ex::grad( up ) && ex::grad( up ) )*mu/rho0; +// (*eps).write("e"); +// +// reuseTmp xit( mesh ); +// std::shared_ptr > xi( xit() ); +// *xi = 2.0*( ex::grad( Tp ) & ex::grad( Tp ) )*ka/rho0; +// (*xi).write("x"); +// +// reuseTmp vortt( mesh ); +// std::shared_ptr > vort( vortt() ); +// *vort = ex::curl( U ); +// (*vort).write("v"); +// +// reuseTmp gpt( mesh ); +// std::shared_ptr > gp( gpt() ); +// *gp = ex::grad(p); +// (*gp).write("z"); +// +// mesh.writeSOS(); +// } + } + + + std::cout<< timer.elapsed().wall / 1e9 <::update for( int j=0; jnj(); j++ ) { int k=settings::m()/2; - fld_ptr->operator()(i,j,k) = this->val(); - } + component(fld_ptr->operator()(i, j, k), this->comp()) = this->val(); + + for( int l=0; loperator()(i, j, k-(l+1)), this->comp()) = 2.0*(l+1)*this->val() - component(fld_ptr->operator()(i, j, k+(l+1) ), this->comp()); + } + + } } } else if( this->dir() == top && parallelCom::k() == parallelCom::nk()-1 ) @@ -22,9 +28,14 @@ void dirichletBC::update for( int j=0; jnj(); j++ ) { int k=fld_ptr->nk()-settings::m()/2-1; - fld_ptr->operator()(i,j,k) = this->val(); + component(fld_ptr->operator()(i, j, k), this->comp()) = this->val(); + + for( int l=0; loperator()(i, j, k+(l+1)), this->comp()) = 2.0*(l+1)*this->val() - component(fld_ptr->operator()(i, j, k-(l+1) ), this->comp()); + } + } } } - } diff --git a/src/BC/neumannBC.hxx b/src/BC/neumannBC.hxx index 947b45f..a9479c7 100644 --- a/src/BC/neumannBC.hxx +++ b/src/BC/neumannBC.hxx @@ -4,27 +4,43 @@ void neumannBC::update Field* fld_ptr ) { - if( this->dir() == bottom && parallelCom::k() == 0 ) - { - for( int i=0; ini(); i++ ) + if( this->dir()== bottom && parallelCom::k() == 0 ) { - for( int j=0; jnj(); j++ ) - { - int k=settings::m()/2; - fld_ptr->operator()(i,j,k) = this->val(); - } - } - } - else if( this->dir() == top && parallelCom::k() == parallelCom::nk()-1 ) - { - for( int i=0; ini(); i++ ) + for( int i=0; ini(); i++ ) + { + for( int j=0; jnj(); j++ ) + { + int k=settings::m()/2; + + component(fld_ptr->operator()(i, j, k), this->comp()) = component(fld_ptr->operator()(i, j, k+1), this->comp()) + -this->val()*fld_ptr->mesh().dz(); + + for( int l=0; loperator()(i, j, k-(l+1)), this->comp()) = component(fld_ptr->operator()(i, j, k+(l+1)), this->comp()) + -2.0*(l+1)*this->val()*fld_ptr->mesh().dz(); + } + } + } + } + + else if( this->dir() == top && parallelCom::k() == parallelCom::nk()-1 ) { - for( int j=0; jnj(); j++ ) - { - int k=fld_ptr->nk()-settings::m()/2-1; - fld_ptr->operator()(i,j,k) = this->val(); - } - } - } -} + for( int i=0; ini(); i++ ) + { + for( int j=0; jnj(); j++ ) + { + int k=fld_ptr->nk()-settings::m()/2-1; + + component(fld_ptr->operator()(i, j, k), this->comp()) = component(fld_ptr->operator()(i, j, k-1), this->comp()) + + this->val()*fld_ptr->mesh().dz(); + for( int l=0; loperator()(i, j, k+(l+1)), this->comp()) = component(fld_ptr->operator()(i, j, k-(l+1)), this->comp()) + -2.0*(l+1)*this->val()*fld_ptr->mesh().dz(); + } + } + } + } +} diff --git a/src/Mesh/Mesh.cpp b/src/Mesh/Mesh.cpp index 4c1dd5e..b6440f8 100644 --- a/src/Mesh/Mesh.cpp +++ b/src/Mesh/Mesh.cpp @@ -103,7 +103,7 @@ Mesh::Mesh( std::string fileName, Time& runTime ) } origin_ = vector( (ni_-1)*dx_*parallelCom::i()+ox, (nj_-1)*dy_*parallelCom::j()+oy, (nk_-1)*dz_*parallelCom::k()+oz ); -#ifdef HAVE_FFTMPI +#ifdef HAVE_HEFFTE localStart_[0] = origin_[0]; localStart_[1] = origin_[1]; localStart_[2] = origin_[2]; @@ -326,6 +326,7 @@ void Mesh::writeCase( std::string path ) fout << "vector per node: velocity " << boost::format("U.%|04|.*****.dat") % parallelCom::myProcNo() << std::endl ; fout << "scalar per node: p " << boost::format("p.%|04|.*****.dat") % parallelCom::myProcNo() << std::endl ; fout << "scalar per node: T " << boost::format("T.%|04|.*****.dat") % parallelCom::myProcNo() << std::endl ; + fout << "vector per node: B " << boost::format("B.%|04|.*****.dat") % parallelCom::myProcNo() << std::endl ; fout << "TIME" << std::endl ; fout << "time set: 1" << std::endl ; @@ -370,7 +371,7 @@ scalar Mesh::z( int k ) const if( hyperbollic_ ) { scalar A = std::pow( s1_/s2_, 0.5 ); -#ifdef HAVE_FFTMPI +#ifdef HAVE_HEFFTE int globK = k-settings::m()/2 + localStart_[2]; #else int globK = k-settings::m()/2 + (parallelCom::k())*(nk_-settings::m()-1); diff --git a/src/Mesh/Mesh.h b/src/Mesh/Mesh.h index 17794a3..944eb6a 100644 --- a/src/Mesh/Mesh.h +++ b/src/Mesh/Mesh.h @@ -67,7 +67,7 @@ class Mesh bool hyperbollic_; -#ifdef HAVE_FFTMPI +#ifdef HAVE_HEFFTE ptrdiff_t localStart_[3]; #endif @@ -156,7 +156,7 @@ class Mesh return bboxHalo_; } -#ifdef HAVE_FFTMPI +#ifdef HAVE_HEFFTE ptrdiff_t* glob_n() { return glob_n_; diff --git a/src/parallelCom/parallelCom.cpp b/src/parallelCom/parallelCom.cpp index 7911b28..4f16883 100644 --- a/src/parallelCom/parallelCom.cpp +++ b/src/parallelCom/parallelCom.cpp @@ -60,13 +60,13 @@ void parallelCom::decompose(int ni, int nj, int nk) { int np=worldSize_; - double lsq=1e10; + double lsq=1e10; for( int ii=1; ii<=np; ii++ ) { for( int jj=1; jj<=np/ii; jj++ ) { - if( std::fabs( double(np)/double(ii)/double(jj) - np/ii/jj ) < 1e-5 ) + if( std::fabs( double(np)/double(ii)/double(jj) - np/ii/jj ) < 1e-5 ) { int kk=np/ii/jj; @@ -82,6 +82,28 @@ void parallelCom::decompose(int ni, int nj, int nk) } } } +// pencil decomposition template + +// int ii=1; +// for( int jj=1; jj<=np; jj++ ) +// { +// if( std::fabs( double(np)/double(jj) - np/jj ) < 1e-5 ) +// { +// int kk=np/jj; +// +// double lsqc = std::pow( double(nj)/double(jj), 2 ) + std::pow( double(nk)/double(kk), 2 ); +// +// if( lsqc < lsq ) +// { +// lsq=lsqc; +// ni_=ii; +// nj_=jj; +// nk_=kk; +// } +// } +// } + + if( master() ) { diff --git a/src/poisson/poisson.cpp b/src/poisson/poisson.cpp index 4acbe98..77c1ccf 100644 --- a/src/poisson/poisson.cpp +++ b/src/poisson/poisson.cpp @@ -12,22 +12,17 @@ poisson::poisson( std::shared_ptr > pptr ) nk_(pptr->mesh().nk()-settings::m()-1), localStart_(pptr->mesh().localStart()) { -fft = std::make_unique(parallelCom::world(), 2); - //fft options - //tflag=0; - permute = 0; - - NFAST = n_[0]; - NMID = n_[1]; - NSLOW = n_[2]; - -// int MPI_Init( int *argc, char ***argv ); - MPI_Comm world = parallelCom::world(); - + + int NFAST = n_[0]; + int NMID = n_[1]; + int NSLOW = n_[2]; + nfast = NFAST; nmid = NMID; nslow = NSLOW; + MPI_Comm world = parallelCom::world(); + npfast = parallelCom::ni(); npmid = parallelCom::nj(); npslow = parallelCom::nk(); @@ -69,9 +64,9 @@ fft = std::make_unique(parallelCom::world(), 2); iglobal = inxlo + ilocal; jglobal = inylo + jlocal; kglobal = inzlo + klocal; - printf("Value (%d,%d,%d) on proc %d \n", - iglobal,jglobal,kglobal, - parallelCom::myProcNo()); +// printf("Value (%d,%d,%d) on proc %d \n", +// iglobal,jglobal,kglobal, +// parallelCom::myProcNo()); if (m == 0) { localStart_[0] = iglobal; @@ -83,28 +78,28 @@ fft = std::make_unique(parallelCom::world(), 2); if (parallelCom::myProcNo() < parallelCom::worldSize()) MPI_Send(&tmp,0,MPI_INT,parallelCom::myProcNo()+1,0,world); } } + //sanity check for fft coordinates + printf("Assigned local start coordinates (%d,%d,%d) on proc %d \n", + localStart_[0],localStart_[1],localStart_[2],parallelCom::myProcNo()); + printf("local range ijk (%d,%d,%d,%d,%d,%d)\n", inxlo,inxhi,inylo,inyhi,inzlo,inzhi); -std::cout<setup(nfast, nmid, nslow, // 3d verion - ilo, ihi, jlo, - jhi, klo, khi, - ilo, ihi, jlo, - jhi, klo, khi, - permute, fftsize, sendsize, recvsize); + //initialize local fft domain + up_N = std::make_unique>(std::array ({inxlo, inylo, inzlo}) ,std::array ({inxhi, inyhi, inzhi})); + N = up_N.get(); - //fft allocate memory - std::cout<(2*fftsize); - phi_ = up_phi_.get(); + //allocate memory for fft + fft = std::make_unique>(*N, *N, world); - up_k2_ = std::make_unique(2*fftsize); + up_k2_ = std::make_unique((fft->size_inbox())); k2_ = up_k2_.get(); + up_phi_ = std::make_unique>((fft->size_inbox())); + phi_ = up_phi_.get(); + + up_phihat_ = std::make_unique>>((fft->size_outbox())); + phihat_ = up_phihat_.get(); + wavenumbers_(); } @@ -118,51 +113,39 @@ void poisson::rhs( std::shared_ptr > f ) { for( int i=(settings::m()/2); i<(pptr_->ni()-1-settings::m()/2); i++ ) { - phi_[l] = (double) f->operator()(i, j, k); - l++; - phi_[l] = 0.0; + phi_[0][l] = f->operator()(i, j, k) ; l++; } - } - } - fft->compute(phi_,phi_,1); + } + } } void poisson::solve() { - int l=0; - for( int k=0; kforward(phi_->data(), phihat_->data()); + + int l=0; + for( int k=0; k(0.0, 0.0); } - else + else { - phi_[l] = -phi_[l] / k2_[l]; + phihat_[0][l] = std::complex(( -std::real(phihat_[0][l]) / k2_[l]) , -std::imag(phihat_[0][l]) / k2_[l]); } l++; - - if( k2_[l] < tools::eps ) - { - phi_[l] = 0.0; - } - else - { - phi_[l] = -phi_[l] / k2_[l]; - } - l++; - } } } - - fft->compute(phi_,phi_,-1); + + fft->backward(phihat_->data(), phi_->data()); l=0; for( int k=(settings::m()/2); k<(pptr_->nk()-1-settings::m()/2); k++ ) @@ -171,14 +154,12 @@ void poisson::solve() { for( int i=(settings::m()/2); i<(pptr_->ni()-1-settings::m()/2); i++ ) { - pptr_->operator()(i,j,k) = phi_[l]; + pptr_->operator()(i,j,k) = phi_[0][l] / (n_[0]*n_[1]*n_[2]); l++; - l++; } } } - } void poisson::wavenumbers_() @@ -186,17 +167,15 @@ void poisson::wavenumbers_() Mesh& m = pptr_->mesh(); auto pi = tools::pi; - printf("Assigned local start coordinates (%d,%d,%d) on proc %d \n", - localStart_[0],localStart_[1],localStart_[2],parallelCom::myProcNo()); scalar kx, ky, kz; //TODO: will be complex for non-central schemes int l=0; - for( int k=0; k -#include "fft3d.h" +#include "heffte.h" + #include "Mesh/Mesh.h" #include "Field/Field.h" #include "Tools/Tools.h" #include #include +#include -using namespace FFTMPI_NS; +using namespace heffte; class poisson { protected: - std::unique_ptr fft; + std::unique_ptr> fft; std::shared_ptr > pptr_; int fftsize; unsigned ni_, nj_, nk_; @@ -51,16 +53,22 @@ class poisson int nfft_in; int nfft_out; int ilo,ihi,jlo,jhi,klo,khi; - int permute,sendsize,recvsize; + int permute,exchange,sendsize,recvsize; int flag,niter,tflag; double tmax; int memorysize; + + std::unique_ptr>> up_phihat_; + std::vector>* phihat_; + + std::unique_ptr> up_phi_; + std::vector* phi_; - std::unique_ptr up_phi_; - FFT_SCALAR* phi_; std::unique_ptr up_k2_; double* k2_; - + + heffte::box3d<>* N; + std::unique_ptr> up_N; public: poisson( std::shared_ptr > ); @@ -74,4 +82,3 @@ class poisson #endif - From 3c6b1799a3afc2dcdf270cc8b413fc4ada77963e Mon Sep 17 00:00:00 2001 From: u43070nw Date: Tue, 28 Jan 2025 16:52:10 +0000 Subject: [PATCH 2/3] Added HDF5 file writing in parallel. --- .DS_Store | Bin 0 -> 6148 bytes CMakeLists.txt | 6 + app/.DS_Store | Bin 0 -> 6148 bytes app/CMakeLists.txt | 15 +- app/OT_Poisson/OT.cpp | 51 ++- app/OT_Poisson/{RKLoop.H => UEqn.H} | 7 +- app/TG/TG.cpp | 2 +- app/TG_Piosson/RKLoop.H | 12 +- app/TG_Piosson/TG.cpp | 3 +- app/mhdjet/BEqn.H | 56 +++ app/mhdjet/UEqn.H | 59 +++ app/mhdjet/mhdjet.cpp | 472 +++++++++++++++++++++ app/shearLayer/BEqn.H | 56 +++ app/shearLayer/UEqn.H | 45 +- app/shearLayer/hydrostaticP.H | 39 -- app/shearLayer/layerAverage.H | 52 --- app/shearLayer/pEqn.H | 4 - app/shearLayer/rhoEqn.H | 38 -- app/shearLayer/shearLayer.cpp | 134 +++--- examples/.DS_Store | Bin 0 -> 10244 bytes examples/OT/.DS_Store | Bin 0 -> 6148 bytes examples/TG/UBC | 36 +- examples/TG/mesh | 2 +- examples/TG/pBC | 12 +- examples/mhdjet/.DS_Store | Bin 0 -> 10244 bytes examples/mhdjet/BBC | 18 + examples/mhdjet/UBC | 18 + examples/mhdjet/mesh | 3 + examples/mhdjet/pBBC | 6 + examples/{shearLayer/pHydBC => mhdjet/pBC} | 4 +- examples/mhdjet/rhoBC | 6 + examples/shearLayer/BBC | 18 + examples/shearLayer/UBC | 36 +- examples/shearLayer/mesh | 5 +- examples/shearLayer/pBBC | 6 + examples/shearLayer/pBC | 8 +- src/BC/dirichletBC.hxx | 88 +++- src/BC/neumannBC.hxx | 94 +++- src/Field/Field.h | 2 +- src/Mesh/Mesh.cpp | 10 +- src/Mesh/Mesh.h | 5 +- src/poisson/poisson.cpp | 270 +++++++++++- src/poisson/poisson.h | 41 +- src/reuseTmp/reuseTmp.hxx | 3 +- 44 files changed, 1429 insertions(+), 313 deletions(-) create mode 100644 .DS_Store create mode 100644 app/.DS_Store rename app/OT_Poisson/{RKLoop.H => UEqn.H} (88%) create mode 100644 app/mhdjet/BEqn.H create mode 100644 app/mhdjet/UEqn.H create mode 100644 app/mhdjet/mhdjet.cpp create mode 100644 app/shearLayer/BEqn.H delete mode 100644 app/shearLayer/hydrostaticP.H delete mode 100644 app/shearLayer/layerAverage.H delete mode 100644 app/shearLayer/pEqn.H delete mode 100644 app/shearLayer/rhoEqn.H create mode 100644 examples/.DS_Store create mode 100644 examples/OT/.DS_Store create mode 100644 examples/mhdjet/.DS_Store create mode 100644 examples/mhdjet/BBC create mode 100644 examples/mhdjet/UBC create mode 100644 examples/mhdjet/mesh create mode 100644 examples/mhdjet/pBBC rename examples/{shearLayer/pHydBC => mhdjet/pBC} (66%) create mode 100644 examples/mhdjet/rhoBC create mode 100644 examples/shearLayer/BBC create mode 100644 examples/shearLayer/pBBC diff --git a/.DS_Store b/.DS_Store new file mode 100644 index 0000000000000000000000000000000000000000..fd28e08626419453760300e8908c97e8b9d4f55a GIT binary patch literal 6148 zcmeHKOH0E*5Z-O8-BN@e6nb3nTCfivh?fxSA26Z^m736^!8BW%)Ci@Jv;HA}iND90 z-3?gUgC`L?1DkJl9=kgqWFL$%?#_lC#u|(<0S%F(vPRIn)>X2>h+K~`VG$d#B*>D< znt}eJ3D+)Q&9BU7WA<(RL8JokeFT#zOwxAylb6bst?jC)in_S>o>b#9A$&X?zxInKa59Hoe&O35OQ}D#zU34Y8DR@o$DEgsEOL3(VWkF-J_Q5 zopcs0IX~^TTk^1Xyjax4-u}Va<={CQ$Lh_H$$@t*I~EIg2W6$M7jGIzDt-iCkypeM z5(C5lF+dD#4g>ZK5cSO|nJOg)h=Ct6fct}lhUi$#4eG4}K6rgbe+3Z*bbL!7N{f!g z+#q;BxJd;xsoXv>xJd`Qv~iBb+@MKkT+a;S*qMvR3)i!QUFvYg9fR~G28e-W2Fj+} z!1MnCewnq8{N)tt5d*})KVyJb`d;6KqU_naQXZbQ650ba6pYJJ0ResK5&#|CN1EDc c{1SDDb1dctaTK)cbU?ZYXhP^i4EzEEU;6${9smFU literal 0 HcmV?d00001 diff --git a/CMakeLists.txt b/CMakeLists.txt index cb11345..8c9586d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -81,6 +81,12 @@ if(PROFILE) set(LIBS ${LIBS} "${GPERFTOOLS_LIBRARIES}") endif() +find_package(HDF5 COMPONENTS C CXX HL REQUIRED) +include_directories( ${HDF5_INCLUDES} ) +include_directories (${HDF5_INCLUDE_DIR}) +set (LIBS ${LIBS} "${HDF5_CXX_LIBRARIES}") + + include_directories(src) add_subdirectory (src) diff --git a/app/.DS_Store b/app/.DS_Store new file mode 100644 index 0000000000000000000000000000000000000000..d3a5c486434b2765aeb4cf95c8447bc1844da532 GIT binary patch literal 6148 zcmeHKQA@)x5WdW*Eu!o}!N-8F12=UJ@uf_C7c2UpGFv;eSQ}|O_b>*1)<5Jg@%MO_ zWa60mB8ZGTxO|t(T|&N;Tmt~2(~p_}H2`o>2`er(UkHtpPD#aj2!);_g&YdVpbs}w z(d_t-4A9xF;EDC{#1wo!zXFaSfGiy3k71%hjC&0qlQ=ILjWB&dP@*Wi=M^OAda92lZt3kg?(ZOla79A<2;L{L6Z(buZ;89m4$tw2)#P`r49$- z8RV83U jLziNVrBb|(ss;U$3`EajX%Iaq{3D=g;D#CaQwH7vIT=x6 literal 0 HcmV?d00001 diff --git a/app/CMakeLists.txt b/app/CMakeLists.txt index 43d9daa..663d8bb 100644 --- a/app/CMakeLists.txt +++ b/app/CMakeLists.txt @@ -1,19 +1,24 @@ add_executable (TG ./TG/TG.cpp) -target_link_libraries (TG NOVA ${Boost_LIBRARIES} ${MPI_LIBRARIES} ${FFTMPI_LIBRARIES} ${LIBS}) +target_link_libraries (TG NOVA ${Boost_LIBRARIES} ${MPI_LIBRARIES} ${HEFFTE_LIBRARIES} ${LIBS}) install(TARGETS TG DESTINATION ${CMAKE_INSTALL_PREFIX}/bin) add_executable (TG_Poisson ./TG_Piosson/TG.cpp) -target_link_libraries (TG_Poisson NOVA ${Boost_LIBRARIES} ${MPI_LIBRARIES} ${FFTMPI_LIBRARIES} ${LIBS}) +target_link_libraries (TG_Poisson NOVA ${Boost_LIBRARIES} ${MPI_LIBRARIES} ${LIBS}) install(TARGETS TG_Poisson DESTINATION ${CMAKE_INSTALL_PREFIX}/bin) add_executable (OT ./OT/OT.cpp) -target_link_libraries (OT NOVA ${Boost_LIBRARIES} ${MPI_LIBRARIES} ${FFTMPI_LIBRARIES} ${LIBS}) +target_link_libraries (OT NOVA ${Boost_LIBRARIES} ${MPI_LIBRARIES} ${HEFFTE_LIBRARIES} ${LIBS}) install(TARGETS OT DESTINATION ${CMAKE_INSTALL_PREFIX}/bin) add_executable (OT_Poisson ./OT_Poisson/OT.cpp) -target_link_libraries (OT_Poisson NOVA ${Boost_LIBRARIES} ${MPI_LIBRARIES} ${FFTMPI_LIBRARIES} ${LIBS}) +target_link_libraries (OT_Poisson NOVA ${Boost_LIBRARIES} ${MPI_LIBRARIES} ${HEFFTE_LIBRARIES} ${LIBS}) install(TARGETS OT DESTINATION ${CMAKE_INSTALL_PREFIX}/bin) add_executable (shearLayer ./shearLayer/shearLayer.cpp) -target_link_libraries (shearLayer NOVA ${Boost_LIBRARIES} ${MPI_LIBRARIES} ${FFTMPI_LIBRARIES} ${LIBS}) +target_link_libraries (shearLayer NOVA ${Boost_LIBRARIES} ${MPI_LIBRARIES} ${LIBS}) install(TARGETS shearLayer DESTINATION ${CMAKE_INSTALL_PREFIX}/bin) + +add_executable (mhdjet ./mhdjet/mhdjet.cpp) +target_link_libraries (mhdjet NOVA ${Boost_LIBRARIES} ${MPI_LIBRARIES} ${LIBS}) +install(TARGETS mhdjet DESTINATION ${CMAKE_INSTALL_PREFIX}/bin) + diff --git a/app/OT_Poisson/OT.cpp b/app/OT_Poisson/OT.cpp index 9943f45..3bd3cd7 100644 --- a/app/OT_Poisson/OT.cpp +++ b/app/OT_Poisson/OT.cpp @@ -23,7 +23,7 @@ int main(int argc, char* argv[]) { settings::process( argc, argv ); - Time time( 0.005, 8.01, 20 ); //args: dt, endT, write interval / steps + Time time( 0.001, 2.01, 10 ); //args: dt, endT, write interval / steps const scalar pi = tools::pi; parallelCom::decompose( settings::zoneName()+"/"+"mesh" ); @@ -45,8 +45,8 @@ int main(int argc, char* argv[]) std::shared_ptr > pB_ptr( std::make_shared >( mesh, 0, "pB" ) ); auto& pB = (*pB_ptr); - scalar mu = 0.01; - scalar eta = 0.01; + scalar mu = 0.057; + scalar eta = 0.057; poisson pEqn(p_ptr); poisson pBEqn(pB_ptr); @@ -54,6 +54,26 @@ int main(int argc, char* argv[]) std::default_random_engine eng( parallelCom::myProcNo() ); std::uniform_real_distribution dist(-1.0,1.0); +/* + //initial conditions + for( int i=settings::m()/2; i() ); reduce( Jmax, maxOp() ); + reduce( Umax, maxOp() ); reduce( Em, plusOp() ); reduce( Jav, plusOp() ); reduce( omega, plusOp() ); @@ -146,15 +170,18 @@ int main(int argc, char* argv[]) Jav /= n; omega /= n; +// Jmax=Jmax*17.54; + scalar epsilon=0.0; epsilon = -mu*(omega*omega) - eta*(Jav*Jav); if( time.curTime() > time.dt() && parallelCom::master() ) { - data< >, 4> dU; for( int i=0; i<4; i++ ) { @@ -25,17 +26,18 @@ { Ustar = Ustar+settings::RK4a()[rk][j]*dU[j]; } - + } Ustar.correctBoundaryConditions(); - pEqn.rhs(ex::div(Ustar) / time.dt() / settings::RK4c()[rk]); + + pEqn.rhs(ex::div(Ustar) / time.dt() ); pEqn.solve(); p.correctBoundaryConditions(); - U = Ustar - settings::RK4c()[rk]*time.dt() * ex::grad(p); + U = Ustar - time.dt() * ex::grad(p); U.correctBoundaryConditions(); - } - + + if( time.writelog() ) { diff --git a/app/TG_Piosson/TG.cpp b/app/TG_Piosson/TG.cpp index 05effcc..8c85fb6 100644 --- a/app/TG_Piosson/TG.cpp +++ b/app/TG_Piosson/TG.cpp @@ -22,7 +22,7 @@ int main(int argc, char* argv[]) { settings::process( argc, argv ); - Time time( 0.01, 20, 10 ); //args: dt, endT, write interval / steps + Time time( 0.01, 20, 1 ); //args: dt, endT, write interval / steps const scalar pi = tools::pi; parallelCom::decompose( settings::zoneName()+"/"+"mesh" ); @@ -53,6 +53,7 @@ int main(int argc, char* argv[]) scalar z = -pi+(k-settings::m()/2)*mesh.dz()+mesh.origin().z(); U(i, j, k) = vector(std::sin(x)*std::cos(y)*std::cos(z),-std::cos(x)*std::sin(y)*std::cos(z),0); +// U(i, j, k) = vector(0.0, 0.0, 0.0); } } } diff --git a/app/mhdjet/BEqn.H b/app/mhdjet/BEqn.H new file mode 100644 index 0000000..f7861ab --- /dev/null +++ b/app/mhdjet/BEqn.H @@ -0,0 +1,56 @@ +{ + std::array >, 4> dB; + for( int i=0; i<4; i++ ) + { + reuseTmp dBt( mesh ); + dB[i]= dBt(); + } + + + reuseTmp Boldt( mesh ); + std::shared_ptr > Bold( Boldt() ); + + *Bold = B; + + for( int rk=0; rk<4; rk++ ) + { + + dB[rk] = + time.dt() * + ( + -ex::div( U, B ) + +ex::div( B, U ) + +ex::laplacian( eta, B ) + ); + + Bstar = Bold; + for( int j=0; j<=rk; j++ ) + { + Bstar = Bstar + settings::RK4a()[rk][j]*dB[j]; + } + + B = Bstar; + B.correctBoundaryConditions(); + + } + + pBEqn.rhs(ex::div(B)); + pBEqn.solve(); + pB.correctBoundaryConditions(); + B -= ex::grad(pB); + B.correctBoundaryConditions(); + + J = ex::curl(B); + J.correctBoundaryConditions(); + + if( time.writelog() ) + { + auto divB = ex::div(B); + auto maxdivB = tools::maximum(divB); + if( parallelCom::master() ) + { + std::cout<<"Magnetic field divergence: " << maxdivB< >, 4> dU; + for( int i=0; i<4; i++ ) + { + reuseTmp dUt( mesh ); + dU[i]= dUt(); + } + + reuseTmp UOldt( mesh ); + std::shared_ptr > UOld( UOldt() ); + + reuseTmp Ustartt( mesh ); + std::shared_ptr > Ustart( Ustartt() ); + *Ustart = Ustart; + + *UOld = U; + + for( int rk=0; rk<4; rk++ ) + { + dU[rk] = + time.dt() * + ( + -ex::div( U, U ) + +ex::laplacian( mu, U ) + +ex::div( B, B ) + -ex::grad( 0.5*(B&B) ) + ); + + Ustar = UOld; + for( int j=0; j<=rk; j++ ) + { + Ustar = Ustar+settings::RK4a()[rk][j]*dU[j]; + } + + } + + Ustar.correctBoundaryConditions(); + pEqn.rhs(ex::div(Ustar) / time.dt()); + pEqn.solve(); + p.correctBoundaryConditions(); + U = Ustar - time.dt() * ex::grad(p); + + U.correctBoundaryConditions(); + + Umag = (U&U); + + if( time.writelog() ) + { + auto divU = ex::div(U); + auto maxdivU = tools::maximum(divU); + + if( parallelCom::master() ) + { + std::cout<<"Velocity divergence: " << maxdivU< +#include "Types/vector/vector.h" +#include "ex/grad/grad.h" +#include "ex/div/div.h" +#include "ex/curl/curl.h" +#include "ex/laplacian/laplacian.h" +#include "Time/Time.h" +#include "parallelCom/parallelCom.h" +#include "settings/settings.h" +#include "BC/BC.h" +#include +#include "Tools/Tools.h" +#include "poisson/poisson.h" +#include + +#include "H5cpp.h" + +#include +#include +#include + + //H5 Create + //H5 Write dataset using MPI + void writeVectorFieldToHDF5( MPI_Comm comm, const std::string &fileName, const std::vector>> &vectorField) { + + MPI_Comm world = parallelCom::world(); + int rank, size; + MPI_Comm_rank(comm, &rank); + MPI_Comm_size(comm, &size); + + // Dimensions of the vector field + hsize_t hni = vectorField.size(); + hsize_t hnj = vectorField[0].size(); + hsize_t hnk = vectorField[0][0].size(); + + //prepare data + std::vector data(hni * hnj * hnk ); + int n=0; + for (int i = 0; i < hni; ++i) { + for (int j = 0; j < hnj; ++j) { + for (int k = 0; k < hnk; ++k) { + data[n] = vectorField[i][j][k]; + n++; + } + } + } + + // Create an HDF5 file + hid_t plist_id = H5Pcreate(H5P_FILE_ACCESS); + H5Pset_fapl_mpio(plist_id, parallelCom::world(), MPI_INFO_NULL); + hid_t file_id = H5Fcreate(fileName.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, plist_id); + + // Create a data space + hsize_t dims[] = { static_cast((hni)*parallelCom::ni() ), + static_cast((hnj)*parallelCom::nj() ), + static_cast((hnk)*parallelCom::nk() ) }; // 3 for the vector components (x, y, z) + hid_t dataspace_id = H5Screate_simple(3, dims, NULL); + + //Create a dataset + hid_t dataset_id = H5Dcreate2(file_id, "/3DArray", H5T_NATIVE_DOUBLE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + + // Create the local dataspace for each MPI rank + hsize_t local_dims[3] = {hni, hnj, hnk}; + hsize_t start[3] = {static_cast((hni*parallelCom::i()) ) , static_cast((hnj*parallelCom::j()) ) , static_cast((hnk*parallelCom::k()) ) }; + hsize_t count[3] = {hni, hnj, hnk}; + hid_t memspace_id = H5Screate_simple(3, local_dims, NULL); + + + //Sanity Checks for data indexing +// std::cout << " ni: " << parallelCom::ni() << " nj: " << parallelCom::nj() << " nk: " << parallelCom::nk() << " i: " < &vectorField) { + // Dimensions of the vector field + hsize_t hni = vectorField[0]-settings::m(); + hsize_t hnj = vectorField[1]-settings::m(); + hsize_t hnk = vectorField[2]-settings::m(); + + // Create an HDF5 file + H5::H5File file(fileName, H5F_ACC_TRUNC); + + // Create a data space + hsize_t dims[3] = {hni, hnj, hnk }; // 3 for the vector components (x, y, z) + H5::DataSpace dataspace(3, dims); + + // Create a dataset + H5::DataSet dataset = file.createDataSet("VectorField", H5::PredType::IEEE_F64LE, dataspace); + + // Prepare data to write + std::vector data(hni * hnj * hnk * 3); + for (int i = 0; i < hni; ++i) { + for (int j = 0; j < hnj; ++j) { + for (int k = 0; k < hnk; ++k) { + data[(i * hnj * hnk + j * hnk + k) + 0] = vectorField[0]; // x component + data[(i * hnj * hnk + j * hnk + k) + 1] = vectorField[1]; // y component + data[(i * hnj * hnk + j * hnk + k) + 2] = vectorField[2]; // z component + + } + } + } + + // Write data to the dataset + dataset.write(data.data(), H5::PredType::IEEE_F64LE); + + // Close the dataset and file + dataset.close(); + file.close(); +} +*/ + + + +int main(int argc, char* argv[]) +{ + settings::restart() = true; + + settings::process( argc, argv ); + Time time( 0.02, 160.02, 50 ); //args: dt, endT, write interval / steps + + const scalar pi = tools::pi; + parallelCom::decompose( settings::zoneName()+"/"+"mesh" ); + + Mesh mesh( settings::zoneName()+"/"+"mesh", time ); + + mesh.write(settings::zoneName()+"/data"); + + + Field U( mesh, vector(0,0,0), "U" ); + Field Ustar( mesh, vector(0,0,0), "U" ); + Field Ustart( mesh, vector(0,0,0), "U" ); + Field B( mesh, vector(0,0,0), "B" ); + Field Bstar( mesh, vector(0,0,0), "B" ); + Field J( mesh, vector(0,0,0), "B" ); + Field omegav( mesh, vector(0,0,0), "U" ); + + std::shared_ptr > Umag_ptr(std::make_shared >(mesh, 0, "U") ); + auto& Umag = (*Umag_ptr); + + std::shared_ptr > p_ptr( std::make_shared >( mesh, 0, "p" ) ); + auto& p = (*p_ptr); + + std::shared_ptr > pB_ptr( std::make_shared >( mesh, 0, "pB" ) ); + auto& pB = (*pB_ptr); + + std::shared_ptr > fftx_ptr( std::make_shared >( mesh, 0, "U" ) ); + auto& fftx = (*fftx_ptr); + + scalar mu = 0.000285714; + scalar eta = 0.000285714; + + poisson fftxDo(fftx_ptr); + poisson pEqn(p_ptr); + poisson pBEqn(pB_ptr); + + std::default_random_engine eng( parallelCom::myProcNo() ); + std::uniform_real_distribution dist(-1.0,1.0); + + //initial conditions + for( int i=settings::m()/2; i>> vectorField(hni, std::vector>(hnj, std::vector(hnk))); + //write Ux + int n=0; + for (int k = 0; k < hnk; ++k) { + for (int j = 0; j < hnj; ++j) { + for (int i = 0; i < hni; ++i) { + vectorField[i][j][k] = U( (i+settings::m()/2), (j+settings::m()/2), (k+settings::m()/2) ).x(); // x component + n++; + } + } + } + std::string H5title = "./h5data/Ux" + std::to_string( (time.curTime()-time.dt()) ) + ".h5"; + writeVectorFieldToHDF5(parallelCom::world(), H5title, vectorField); + + //write Uy + n=0; + for (int k = 0; k < hnk; ++k) { + for (int j = 0; j < hnj; ++j) { + for (int i = 0; i < hni; ++i) { + vectorField[i][j][k] = U( (i+settings::m()/2), (j+settings::m()/2), (k+settings::m()/2) ).y(); // x component + n++; + } + } + } + H5title = "./h5data/Uy_" + std::to_string( (time.curTime()-time.dt()) ) + ".h5"; + writeVectorFieldToHDF5(parallelCom::world(), H5title, vectorField); + + //write Uz + n=0; + for (int k = 0; k < hnk; ++k) { + for (int j = 0; j < hnj; ++j) { + for (int i = 0; i < hni; ++i) { + vectorField[i][j][k] = U( (i+settings::m()/2), (j+settings::m()/2), (k+settings::m()/2) ).z(); // x component + n++; + } + } + } + H5title = "./h5data/Uz_" + std::to_string( (time.curTime()-time.dt()) ) + ".h5"; + writeVectorFieldToHDF5(parallelCom::world(), H5title, vectorField); + + //write p + n=0; + for (int k = 0; k < hnk; ++k) { + for (int j = 0; j < hnj; ++j) { + for (int i = 0; i < hni; ++i) { + vectorField[i][j][k] = p( (i+settings::m()/2), (j+settings::m()/2), (k+settings::m()/2) ); // x component + n++; + } + } + } + H5title = "./h5data/p_" + std::to_string( (time.curTime()-time.dt()) ) + ".h5"; + writeVectorFieldToHDF5(parallelCom::world(), H5title, vectorField); + + //write Bx + n=0; + for (int k = 0; k < hnk; ++k) { + for (int j = 0; j < hnj; ++j) { + for (int i = 0; i < hni; ++i) { + vectorField[i][j][k] = B( (i+settings::m()/2), (j+settings::m()/2), (k+settings::m()/2) ).x(); // x component + n++; + } + } + } + H5title = "./h5data/Bx_" + std::to_string( (time.curTime()-time.dt()) ) + ".h5"; + writeVectorFieldToHDF5(parallelCom::world(), H5title, vectorField); + + //write By + n=0; + for (int k = 0; k < hnk; ++k) { + for (int j = 0; j < hnj; ++j) { + for (int i = 0; i < hni; ++i) { + vectorField[i][j][k] = B( (i+settings::m()/2), (j+settings::m()/2), (k+settings::m()/2) ).y(); // x component + n++; + } + } + } + H5title = "./h5data/By_" + std::to_string( (time.curTime()-time.dt()) ) + ".h5"; + writeVectorFieldToHDF5(parallelCom::world(), H5title, vectorField); + + //write Bz + n=0; + for (int k = 0; k < hnk; ++k) { + for (int j = 0; j < hnj; ++j) { + for (int i = 0; i < hni; ++i) { + vectorField[i][j][k] = B( (i+settings::m()/2), (j+settings::m()/2), (k+settings::m()/2) ).z(); // x component + n++; + } + } + } + H5title = "./h5data/Bz_" + std::to_string( (time.curTime()-time.dt()) ) + ".h5"; + writeVectorFieldToHDF5(parallelCom::world(), H5title, vectorField); + + //write Jx + n=0; + for (int k = 0; k < hnk; ++k) { + for (int j = 0; j < hnj; ++j) { + for (int i = 0; i < hni; ++i) { + vectorField[i][j][k] = J( (i+settings::m()/2), (j+settings::m()/2), (k+settings::m()/2) ).x(); // x component + n++; + } + } + } + H5title = "./h5data/Jx_" + std::to_string( (time.curTime()-time.dt()) ) + ".h5"; + writeVectorFieldToHDF5(parallelCom::world(), H5title, vectorField); + + //write Jy + n=0; + for (int k = 0; k < hnk; ++k) { + for (int j = 0; j < hnj; ++j) { + for (int i = 0; i < hni; ++i) { + vectorField[i][j][k] = J( (i+settings::m()/2), (j+settings::m()/2), (k+settings::m()/2) ).y(); // x component + n++; + } + } + } + H5title = "./h5data/Jy_" + std::to_string( (time.curTime()-time.dt()) ) + ".h5"; + writeVectorFieldToHDF5(parallelCom::world(), H5title, vectorField); + + //write Jz + n=0; + for (int k = 0; k < hnk; ++k) { + for (int j = 0; j < hnj; ++j) { + for (int i = 0; i < hni; ++i) { + vectorField[i][j][k] = J( (i+settings::m()/2), (j+settings::m()/2), (k+settings::m()/2) ).z(); // x component + n++; + } + } + } + H5title = "./h5data/Jz_" + std::to_string( (time.curTime()-time.dt()) ) + ".h5"; + writeVectorFieldToHDF5(parallelCom::world(), H5title, vectorField); + + + + + + if( parallelCom::master() ) + { + std::cout << "Vector field written to vector_field.h5" << std::endl; + } + } + + + scalar Ek=0.0; + scalar Em=0.0; + scalar Jmax=0.0; + scalar Umax=0.0; + omegav = ex::curl(U); + scalar Jav=0.0; + scalar omega=0.0; + int n=0; + + for( int i=settings::m()/2; i() ); + reduce( Jmax, maxOp() ); + reduce( Umax, maxOp() ); + reduce( Em, plusOp() ); + reduce( Jav, plusOp() ); + reduce( omega, plusOp() ); + reduce( n, plusOp() ); + + Ek /= n; + Em /= n; + Jav /= n; + omega /= n; + + scalar epsilon=0.0; + epsilon = -mu*(omega*omega) - eta*(Jav*Jav); + + if( time.curTime() > time.dt() && parallelCom::master() ) + { + data< >, 4> dB; + for( int i=0; i<4; i++ ) + { + reuseTmp dBt( mesh ); + dB[i]= dBt(); + } + + + reuseTmp Boldt( mesh ); + std::shared_ptr > Bold( Boldt() ); + + *Bold = B; + + for( int rk=0; rk<4; rk++ ) + { + + dB[rk] = + time.dt() * + ( + -ex::div( U, B ) + +ex::div( B, U ) + +ex::laplacian( eta, B ) + ); + + Bstar = Bold; + for( int j=0; j<=rk; j++ ) + { + Bstar = Bstar + settings::RK4a()[rk][j]*dB[j]; + } + + B = Bstar; + B.correctBoundaryConditions(); + + } + + pBEqn.rhs(ex::div(B)); + pBEqn.solve(); + pB.correctBoundaryConditions(); + B -= ex::grad(pB); + B.correctBoundaryConditions(); + +// J = ex::curl(B); +// J.correctBoundaryConditions(); + + if( time.writelog() ) + { + auto divB = ex::div(B); + auto maxdivB = tools::maximum(divB); + if( parallelCom::master() ) + { + std::cout<<"Magnetic field divergence: " << maxdivB< >, 4> dU; -// for( int i=0; i<4; i++ ) -// { - reuseTmp dUt( mesh ); - std::shared_ptr > dU( dUt() ); -// dU[i]= dUt(); -// } + reuseTmp dUt( mesh ); + std::shared_ptr > dU( dUt() ); reuseTmp UOldt( mesh ); std::shared_ptr > UOld( UOldt() ); + reuseTmp Uct( mesh ); + std::shared_ptr > Uc( Uct() ); + *UOld = U; -// std::shared_ptr > gradU( fdc::grad( U ) ); -// gradU->correctBoundaryConditions(); for( int rk=0; rk<4; rk++ ) { @@ -23,29 +19,34 @@ -ex::div( U, U ) +ex::laplacian( mu, U ) // +mu*ex::div( ( dev2(transpose(gradU) ) ) ) - -( (rho0*T)*g ) +// -( (rho0*T)*(g) ) + +ex::div( B, B ) + -ex::grad( 0.5*(B&B) ) ); + Ustar = UOld; -// for( int j=0; j<=rk; j++ ) -// { -// Ustar = Ustar+settings::RK4[rk]*dU; - Ustar = Ustar+RK4[rk]*dU; -// } + Ustar = Ustar+RK4[rk]*dU; Ustar.correctBoundaryConditions(); -// pEqn.rhs(ex::div(Ustar) / time.dt() / settings::RK4c()[rk]); -// pEqn.solve(); + pEqn.rhs(ex::div(Ustar) / time.dt() / settings::RK4c()[rk]); + pEqn.solve(); p.correctBoundaryConditions(); -// U = Ustar - settings::RK4c()[rk]*time.dt() * ex::grad(p); //(p_rgh?) + U = Ustar - settings::RK4c()[rk]*time.dt() * ex::grad(p); U.correctBoundaryConditions(); } - - - - +// if( time.writelog() ) +// { + auto divU = ex::div(U); + auto maxdivU = tools::maximum(divU); + + if( parallelCom::master() ) + { + std::cout<<"Velocity divergence: " << maxdivU< cumTot( parallelCom::worldSize(), scalar(0.0) ); - cumTot[parallelCom::myProcNo()] = pHydc; - - all_reduce( cumTot, std::plus() ); - - pHydc = 0.0; - double aboveProcs=0.0; - - for( int pk=parallelCom::k()+1; pk=settings::m()/2; k-- ) - { - pHydc += T(i, j, k)*mesh.dz()*g[2]*rho0; - pHyd(i, j, k) = pHydc; // + aboveProcs; - } - - } - } -//add field and implement null boundary condition -// pHyd.correctBoundaryConditions(); - - p += pHyd; -} - diff --git a/app/shearLayer/layerAverage.H b/app/shearLayer/layerAverage.H deleted file mode 100644 index 0e7087e..0000000 --- a/app/shearLayer/layerAverage.H +++ /dev/null @@ -1,52 +0,0 @@ -{ - for( int k=0; k proc_avgU( parallelCom::worldSize(), 0.0 ); - std::vector proc_avgT( parallelCom::worldSize(), 0.0 ); - std::vector proc_N( parallelCom::worldSize(), 0 ); - - for( int j=settings::m()/2; j()); - all_reduce(proc_avgT, std::plus()); - all_reduce(proc_N, std::plus()); - - vector avgU(0,0,0); - scalar avgT=0.0; - int N=0; - - for( int proci=0; procicorrectBoundaryConditions(); - //Tp->correctBoundaryConditions(); - } -} - diff --git a/app/shearLayer/pEqn.H b/app/shearLayer/pEqn.H deleted file mode 100644 index 63303fe..0000000 --- a/app/shearLayer/pEqn.H +++ /dev/null @@ -1,4 +0,0 @@ -{ - p = (c*c/2.0)*rho0*rho0; -} - diff --git a/app/shearLayer/rhoEqn.H b/app/shearLayer/rhoEqn.H deleted file mode 100644 index 4d4334c..0000000 --- a/app/shearLayer/rhoEqn.H +++ /dev/null @@ -1,38 +0,0 @@ -{ - reuseTmp drhot( mesh ); - std::shared_ptr > drho( drhot() ); - - reuseTmp rhoct( mesh ); - std::shared_ptr > rhoc( rhoct() ); - - reuseTmp rhoOldt( mesh ); - std::shared_ptr > rhoOld( rhoOldt() ); - - *rhoOld = rho; - *rhoc = rho; - - for( int rk=0; rk<4; rk++ ) - { - *drho = - time.dt() * - ( - -ex::div( rho*U ) - ); - - //drho->correctBoundaryConditions(); - - *rhoc = rhoc+RK4[rk]*drho; - - if( rk<3 ) - { - rho = rhoOld + (int(rk/2)+1)/2.0 * drho; - } - else - { - rho=rhoc; - } - - rho.correctBoundaryConditions(); - } - -} diff --git a/app/shearLayer/shearLayer.cpp b/app/shearLayer/shearLayer.cpp index 6950223..4b34c15 100644 --- a/app/shearLayer/shearLayer.cpp +++ b/app/shearLayer/shearLayer.cpp @@ -3,7 +3,6 @@ #include "Types/vector/vector.h" #include "ex/grad/grad.h" #include "ex/div/div.h" -#include "ex/curl/curl.h" #include "ex/laplacian/laplacian.h" #include "Time/Time.h" #include "parallelCom/parallelCom.h" @@ -22,14 +21,15 @@ int main(int argc, char* argv[]) { settings::process( argc, argv ); - Time time( 0.1, 1000000, 100 ); //args: dt, endT, write interval / steps + Time time( 0.001, 420, 1 ); //args: dt, endT, write interval / steps + const scalar pi = tools::pi; parallelCom::decompose( settings::zoneName()+"/"+"mesh" ); Mesh mesh( settings::zoneName()+"/"+"mesh", time ); mesh.write(settings::zoneName()+"/data"); - +/*-------------------- scalar Mach=0.2; scalar Re0=1967.0; scalar Ri0=0.08; @@ -43,70 +43,102 @@ int main(int argc, char* argv[]) scalar t0=0.1;//-Ri0*u0*u0/g[2]/h0; scalar mu = rho0*u0*h0/Re0; scalar ka = mu/Pr; - vector g( 0, 0, -Ri0*u0*u0/t0/h0 ); + scalar eta = mu; +// vector g( 0, 0, -Ri0*u0*u0/t0/h0 ); + vector g( 0, 0, -0.000011129 ); Field U( mesh, vector(1,0,0), "U" ); Field T( mesh, 0, "T" ); Field pHyd( mesh,0, "pHyd" ); Field Ustar( mesh, vector(0,0,0), "U" ); -//------------------------------------------- -// time.dt() = 0.5*std::min(mesh.dx(), std::min( mesh.dy(), mesh.dz() ) )/(u0/Mach); -// -// if( parallelCom::master() ) -// { -// std::cout<<"Speed of sound: "< B( mesh, vector(0,0,0), "B" ); +// Field Bstar( mesh, vector(0,0,0), "B" ); + +----------------------*/ + + Field U( mesh, vector(0,0,0), "U" ); + Field Ustar( mesh, vector(0,0,0), "U" ); + Field B( mesh, vector(0,0,0), "B" ); + Field Bstar( mesh, vector(0,0,0), "B" ); + + scalar mu = 1.0/28000.0; + scalar u0=1.0; + scalar delta = 1.0/28.0; + scalar cn = 0.001; + scalar eta = mu*10.0; - scalar drhodz=(std::pow( ( (c*c/2.0/rho0)*rho0*rho0 + (-t0/2.0*g[2]*rho0)*0.1e-8 ) / (c*c/2.0/rho0), 0.5 ) - rho0)/0.1e-8; std::default_random_engine eng( parallelCom::myProcNo() ); std::uniform_real_distribution dist(-1.0,1.0); - std::shared_ptr > p_ptr( std::make_shared >( mesh, 0, "p" ) ); + std::shared_ptr > p_ptr( std::make_shared >( mesh, (0.87), "p" ) ); auto& p = (*p_ptr); + std::shared_ptr > pB_ptr( std::make_shared >( mesh, 0, "pB" ) ); + auto& pB = (*pB_ptr); + poisson pEqn(p_ptr); + poisson pBEqn(pB_ptr); { - #include "pEqn.H" //initial conditions for( int i=settings::m()/2; i6~5m%ac0Og9=lGPY_r}?HjPcPX?E?z#Bu6w{Kx~_iOH^CiCbqMlX#l- zj=QsKH#MOy2-MP6JPHU9P!-BcK*TH13er*}s)|yfq7bP;l&VOON`MeVMdeXA_ukpW zUfVxPA8MpoY3`Z(ICtipIp6v2y~`LwdodGbEXWvB>EhDRsaW90`}sBRr?@Q|NsvBc zCM%?-oQYJ{PR$$bHFzHIJm7i2^ML08&jW7<4-jYbpSoO1uk<|NdBF3)B@giY5TuLC zm@8+b6sdy>u>~Nxgqp>V?r9$&bb_voxpGEIp-OYg-2+@zxLXWJb?OiD;Ur_OoRLzh zGe~s?_sMW~C=gF4y`UJ*5R=j?Jr8&unC}6;cdumQY>Z{tt!4Z7M^h77%PB-6=b%(9 zSz5VFU8XKqZ;p@IQ}IHwP;}ap`NO=HuC2mBH>_kTl8c&I310&M8m=6NYlxanz}NyrS;(G z@$7iczE|1-Hy5z0%IvzY+^CF`&Zg)kxxyT4v(99PJgul`D_MKKWUcm&<@!~=fx)5S zL;9#uGf8G0m?&6w&KWTC4#`oY$|3XjrR`i_-sJs*X**jS&pD$;bvl_%7qiKNxhtEE zTX&e$R^_Dgb~f8*I~Jd2c&K3BUf_N6x`qxpcD~RdCL6M7emTl0I_5!>CbcHbPMX@O zEwdcjsTJ3*S$qA)En8dKI`;PTovkveYxG+En!!oSu~J#Hf6B@f#)gwl+RSCF+_6LC z|6IX3YMFVh8lzewmmD|G))>`Qt5*AbL&cPJThS_TulK!EGm>h;Aw%U9Vf_ZpaJ0m63I?ef)tgi!qxEN#DKjgZ zn)U6f(X9GE}}s`!G`&Xu|{R>95-(|hZDzF|Juu1vLFpoU^i@=Jtb8{M1Imk|yP#Uiq<5zqNqBkN>6>8zK7}vftN0qefp6k89>e2!0zbvGIFH}pw|E{e;cxgmE-IBum9kp# zDc33MmFtyyr9o*`b}C)Ueq~S@QB3*zN(pi!CG*_&qa@3b<(>m7-&+`2fxY|Q-Q9EL zAnUmm97|W%uDs^jAP3m?8{031p$yzotj@-O2>Ur;5n(T)qC8#+ao4T!>w$j4g?W29 zQcDOnmT?5AXZh&0>-2!Jl#nB<8|tZANvM(4P2oVmPzbED8i@vUqk^EQ3EgdBs=}?b z?9g>%34v1fwcSKj#`A%!cXZNB9AfGl73$Bi7uf~&SN3lT^%b}pt5A;!nz54tdq4UR zM*{s6+((c?S_HUF;hl$r_u_pN;&)My-;4Y45quOM!)Nh1d|m|l*F~Uz3qQak_z`}L zGk6kD;pcc7zrZs%hu`4>{*0H)!g>GPaK=B%!+A2E&t>i0F>2=!P}l*Mn03_sNGZhM zGOAoPJ-xN1B^+#OY~fVd(nz$qwW+y`DL*XjT9hd{lL}jA5n;S!%UmM#3mzPq#YS%h z`BAbh=ez}}_;P|J0t~Eg*xdN8ow2Tcz3&mpUe{M%BEZYkCjmB_0Eq#^E&)m``WFM? z0>tQF7t&PY^c<{K7tL|0aU5k(@=0dL*h5GhQY}zrL+m)r`z6V)W}fTB)hK zlx*Z9^p!WMWvuA;cXX;medSq(XEEyAtIp#@YVBE=6UiI<*ZKW(IMKIib&Wsn?_ZP? zA0KOt7hhm(5l+5H0FGZ+n={XTx1@1{>;iE*qA^_A(oC(7J4> zhlmXcE-oAD4a9~76_*VIn~4nxp)MN+8i@@Fe=ZvaHqT;1y_MLIuTce4w&``mhAtJh z2Rexg`I1r%y|EGgEIZF$pa{RnUJ)xlgVucwwBFlJ5xpCGu@61y#Q=sdELMIdj$sV9 zQj8a9-8Y33xP#VxAEXGsTde**Ay$76;vtIjhw){K@^9li_%6}n8CnPaM63g!##veg zzKGwKMR0E^cBRC{QUsR~;oM}HtT;B literal 0 HcmV?d00001 diff --git a/examples/OT/.DS_Store b/examples/OT/.DS_Store new file mode 100644 index 0000000000000000000000000000000000000000..5e783bb5b431d4e35df70badcfbef8827e9725a0 GIT binary patch literal 6148 zcmeHKyH3ME5S#-Sk!V6kdB4CPoTBgp@B@UHF2WK?dUyG4W*28__o31aTjCh# zpAH5e0f-%GH$M9;L99j)TjCf6+UO|#53bKxi*gD?6CDxH;wDqjqTXZ!)da9CQh2>v3Be@DcjwhNiuXj zq zWwL44%`WR18o~&K5eOp?Mj(tp7=arl0>s(;P^+a3uVDnj2!s(>i~z3>F*>=72XaA5 zQ99@%wg4n6$t(`KXZ(g_ERgX)E=Wl&bf9!4DP1wxVxV*co z40b3eo=)S(M0dcrl;JguKp27L5#W3GN-!Y>KBS8G@3jRZ903b5;DQ4lb;?1WG@oo{ zarc#r-RI$KaYGfHBp-|Z-bbhZ7EQB>oRf8(Q}jkW$ydV+J(NH45DnncLoL!gu<1ok zrRThvbk$SDRPH#}#^}pLx^SOL zAGNKyBE9Y19|R!%tBrbS&&OxTX;&JDPE$Ml46LN=(HswG+DxpUo}9c{6$ZClQ_ z#m;T-qWkVGz30zs>e}YbeMge_W~UwZA$m7DFeq5SBVBwvOv+o%G3(6^$~ji7FU;}P z!W=cN+894NF_~m6?^#EzY|6K1t)Y}}c}cd$qv?)jT<1{Enx$z@WZZ0i+M$-(Oe&kn zQ)jC`n@!kfEpO7d=6!xQ<+(Y3P;{KM$(P%btd9F1obhef@s3(KkEEtuShY5K^R1gY zdUos@+_V4S<(k!Nm{zCNPt4k$oz7b0bGGT9K9=$_mSft^sVRz{Oa8 zrmf4YVO`_;h+&lWE{rmnPupkmwtpXMR#o|oCgVLC`=)xrqzN2wO+KP2T8t)RgQ_0q z_Q9M*S83K-j8=`M)G?1HBZp|$Hfzk&#_pw`Sh{v;-6}J+@oXw>WuR-!v_gCf(Sujv0rIgvx$OKFoGzZGtrY1`^q7bG%hT zYDq+z;~mmkUyL`p=Onux+F%!qz!W8v19R{QJPMD&6YyoY0MEep;79NxyaX@9&)`?^ zI{XpdfPcc9NH9hntFZ=Gqk*^J2D}x!un%|QF5HbncmVIhLwFdEViNDcG-lAo9D10? zS-c-VgdfI_;K%XP_!)c>KZ_Ue+jtSbgU{nj_%i+&U&G(y8~9iJ8@{Q4q9_$ggR)ME zDeX$HoPknEqGL;Er0gO6T1gZX3dJv%MA6}qJ$v7gB)YPSr%i3$yXx28)ZEhE)!Vmw z5nK|B1^5-iS|EAJuPu-)GO84XGQYIwI6>W-zgn8HDxQ++W!}3frbXEbN(QOkvVqi< zln7G2tu-2D3gw?vZ{Hl%SOsO5CVaQFl8S+{b(^NKO3E|Yw$mUL0wvy0!|<$8uOEJ& zg0I5Y;UYXw__bZNXODgl*V^TX8$?AOsKL2=2vE z+=mGq#|fOoV}xH5Ej)=Xo*@kTgyQ>g0Y88b;KTSRK8_#7j|nt?1)na#=jc*=;t#Hc z&l1wp0)m@|O~KvScbe*YI=YGKdOP1H)s>l95y3^4=cZn5Oib!iWhv;eFh)eSQ9!jFdjlfdiNgVr?I*xWBAcRRm`OBtDp-v9v>1Co5!nOy zAbya5`#6!wAQ>fi|5vY)_6yx9jhp+zs z{~Pbc!#4^e5Jup-MF1;@5<>$tTw#}8I2#pifZ_~-bR rz!rRs5bg26wjiZ&CqLm%Qhv?::update Field* fld_ptr ) { + //bottom if( this->dir() == bottom && parallelCom::k() == 0 ) { for( int i=0; ini(); i++ ) @@ -21,7 +22,10 @@ void dirichletBC::update } } } - else if( this->dir() == top && parallelCom::k() == parallelCom::nk()-1 ) + + + //top + if( this->dir() == top && parallelCom::k() == parallelCom::nk()-1 ) { for( int i=0; ini(); i++ ) { @@ -38,4 +42,86 @@ void dirichletBC::update } } } + + + //east + if( this->dir() == east && parallelCom::i() == 0 ) + { + for( int j=0; jnj(); j++ ) + { + for( int k=0; knk(); k++ ) + { + int i=settings::m()/2; + component(fld_ptr->operator()(i, j, k), this->comp()) = this->val(); + + for( int l=0; loperator()(i-(l+1), j, k), this->comp()) = 2.0*(l+1)*this->val() - component(fld_ptr->operator()(i+(l+1), j, k ), this->comp()); + } + + } + } + } + + + //west + if( this->dir() == west && parallelCom::i() == parallelCom::ni()-1 ) + { + for( int j=0; jnj(); j++ ) + { + for( int k=0; knk(); k++ ) + { + int i=fld_ptr->ni()-settings::m()/2-1; + component(fld_ptr->operator()(i, j, k), this->comp()) = this->val(); + + for( int l=0; loperator()(i+(l+1), j, k), this->comp()) = 2.0*(l+1)*this->val() - component(fld_ptr->operator()(i-(l+1), j, k ), this->comp()); + } + + } + } + } + + + //south + if( this->dir() == south && parallelCom::j() == 0 ) + { + for( int i=0; ini(); i++ ) + { + for( int k=0; knk(); k++ ) + { + int j=settings::m()/2; + component(fld_ptr->operator()(i, j, k), this->comp()) = this->val(); + + for( int l=0; loperator()(i, j-(l+1), k), this->comp()) = 2.0*(l+1)*this->val() - component(fld_ptr->operator()(i, j+(l+1), k ), this->comp()); + } + + } + } + } + + + //north + if( this->dir() == north && parallelCom::j() == parallelCom::nj()-1 ) + { + for( int i=0; ini(); i++ ) + { + for( int k=0; knk(); k++ ) + { + int j=fld_ptr->nj()-settings::m()/2-1; + component(fld_ptr->operator()(i, j, k), this->comp()) = this->val(); + + for( int l=0; loperator()(i, j+(l+1), k), this->comp()) = 2.0*(l+1)*this->val() - component(fld_ptr->operator()(i, j-(l+1), k ), this->comp()); + } + + } + } + } + + } diff --git a/src/BC/neumannBC.hxx b/src/BC/neumannBC.hxx index a9479c7..43386ed 100644 --- a/src/BC/neumannBC.hxx +++ b/src/BC/neumannBC.hxx @@ -4,7 +4,8 @@ void neumannBC::update Field* fld_ptr ) { - if( this->dir()== bottom && parallelCom::k() == 0 ) + //bottom + if( this->dir() == bottom && parallelCom::k() == 0 ) { for( int i=0; ini(); i++ ) { @@ -24,7 +25,8 @@ void neumannBC::update } } - else if( this->dir() == top && parallelCom::k() == parallelCom::nk()-1 ) + //top + if( this->dir() == top && parallelCom::k() == parallelCom::nk()-1 ) { for( int i=0; ini(); i++ ) { @@ -39,8 +41,94 @@ void neumannBC::update { component(fld_ptr->operator()(i, j, k+(l+1)), this->comp()) = component(fld_ptr->operator()(i, j, k-(l+1)), this->comp()) -2.0*(l+1)*this->val()*fld_ptr->mesh().dz(); - } + } } } } + + + //east + if( this->dir() == east && parallelCom::i() == 0 ) + { + for( int j=0; jnj(); j++ ) + { + for( int k=0; knk(); k++ ) + { + int i=settings::m()/2; + + component(fld_ptr->operator()(i, j, k), this->comp()) = component(fld_ptr->operator()(i+1, j, k), this->comp()) + -this->val()*fld_ptr->mesh().dx(); + + for( int l=0; loperator()(i-(l+1), j, k), this->comp()) = component(fld_ptr->operator()(i+(l+1), j, k), this->comp()) + -2.0*(l+1)*this->val()*fld_ptr->mesh().dx(); + } + } + } + } + + //west + if( this->dir() == west && parallelCom::i() == parallelCom::ni()-1 ) + { + for( int j=0; jnj(); j++ ) + { + for( int k=0; knk(); k++ ) + { + int i=fld_ptr->ni()-settings::m()/2-1; + + component(fld_ptr->operator()(i, j, k), this->comp()) = component(fld_ptr->operator()(i-1, j, k), this->comp()) + + this->val()*fld_ptr->mesh().dx(); + + for( int l=0; loperator()(i+(l+1), j, k), this->comp()) = component(fld_ptr->operator()(i-(l+1), j, k), this->comp()) + -2.0*(l+1)*this->val()*fld_ptr->mesh().dx(); + } + } + } + } + + //south + if( this->dir() == south && parallelCom::j() == 0 ) + { + for( int i=0; ini(); i++ ) + { + for( int k=0; knk(); k++ ) + { + int j=settings::m()/2; + + component(fld_ptr->operator()(i, j, k), this->comp()) = component(fld_ptr->operator()(i, j+1, k), this->comp()) + -this->val()*fld_ptr->mesh().dy(); + + for( int l=0; loperator()(i, j-(l+1), k), this->comp()) = component(fld_ptr->operator()(i, j+(l+1), k), this->comp()) + -2.0*(l+1)*this->val()*fld_ptr->mesh().dy(); + } + } + } + } + + //north + if( this->dir() == north && parallelCom::j() == parallelCom::nj()-1 ) + { + for( int i=0; ini(); i++ ) + { + for( int k=0; knk(); k++ ) + { + int j=fld_ptr->nj()-settings::m()/2-1; + + component(fld_ptr->operator()(i, j, k), this->comp()) = component(fld_ptr->operator()(i, j-1, k), this->comp()) + + this->val()*fld_ptr->mesh().dy(); + + for( int l=0; loperator()(i, j+(l+1), k), this->comp()) = component(fld_ptr->operator()(i, j-(l+1), k), this->comp()) + -2.0*(l+1)*this->val()*fld_ptr->mesh().dy(); + } + } + } + } + } diff --git a/src/Field/Field.h b/src/Field/Field.h index 57edb6c..223b8b8 100644 --- a/src/Field/Field.h +++ b/src/Field/Field.h @@ -18,7 +18,6 @@ */ #ifndef FIELD_H #define FIELD_H - #if defined(_OPENMP) #include #endif @@ -129,6 +128,7 @@ class Field void correctBoundaryConditions(); void write(std::string, std::string); +// void H5write(std::string, std::string); void read(std::string, std::string); void read(std::string, std::string, int); void setBC( std::string ); diff --git a/src/Mesh/Mesh.cpp b/src/Mesh/Mesh.cpp index b6440f8..043deca 100644 --- a/src/Mesh/Mesh.cpp +++ b/src/Mesh/Mesh.cpp @@ -18,7 +18,8 @@ */ #include "Mesh/Mesh.h" #include - +#include +#include "H5Cpp.h" Mesh::Mesh( std::string fileName, Time& runTime ) : @@ -324,9 +325,12 @@ void Mesh::writeCase( std::string path ) fout << "VARIABLE" << std::endl; fout << "vector per node: velocity " << boost::format("U.%|04|.*****.dat") % parallelCom::myProcNo() << std::endl ; + fout << "vector per node: Uprime " << boost::format("Uprime.%|04|.*****.dat") % parallelCom::myProcNo() << std::endl ; + fout << "vector per node: Udash " << boost::format("Udash.%|04|.*****.dat") % parallelCom::myProcNo() << std::endl ; fout << "scalar per node: p " << boost::format("p.%|04|.*****.dat") % parallelCom::myProcNo() << std::endl ; fout << "scalar per node: T " << boost::format("T.%|04|.*****.dat") % parallelCom::myProcNo() << std::endl ; fout << "vector per node: B " << boost::format("B.%|04|.*****.dat") % parallelCom::myProcNo() << std::endl ; + fout << "vector per node: J " << boost::format("J.%|04|.*****.dat") % parallelCom::myProcNo() << std::endl ; fout << "TIME" << std::endl ; fout << "time set: 1" << std::endl ; @@ -344,7 +348,7 @@ void Mesh::writeCase( std::string path ) fout.close() ; - this->writeRestartInfo( path ); + this->writeRestartInfo( path ); } void Mesh::writeRestartInfo( std::string path ) @@ -365,6 +369,7 @@ void Mesh::writeRestartInfo( std::string path ) } + scalar Mesh::z( int k ) const { scalar z; @@ -439,3 +444,4 @@ void Mesh::sendMUIPoints() } } #endif + diff --git a/src/Mesh/Mesh.h b/src/Mesh/Mesh.h index 944eb6a..c655d10 100644 --- a/src/Mesh/Mesh.h +++ b/src/Mesh/Mesh.h @@ -28,6 +28,7 @@ #include #include #include +#include "H5Cpp.h" class Time; @@ -57,7 +58,7 @@ class Mesh vector origin_; - Time& runTime_; +// Time& runTime_; Box bbox_; Box bboxHalo_; @@ -74,6 +75,7 @@ class Mesh public: Mesh( int, int, int, scalar, scalar, scalar, Time& ); Mesh( std::string, Time& ); + Time& runTime_; //access int ni() const @@ -172,6 +174,7 @@ class Mesh void writeCase( std::string ); void writeSOS(); void writeRestartInfo( std::string ); +// void H5write( std::string ); }; #endif diff --git a/src/poisson/poisson.cpp b/src/poisson/poisson.cpp index 77c1ccf..4833fa7 100644 --- a/src/poisson/poisson.cpp +++ b/src/poisson/poisson.cpp @@ -12,7 +12,7 @@ poisson::poisson( std::shared_ptr > pptr ) nk_(pptr->mesh().nk()-settings::m()-1), localStart_(pptr->mesh().localStart()) { - + //Define FFT grid. This can be trimmed now HeFFTe is the prefered library. int NFAST = n_[0]; int NMID = n_[1]; int NSLOW = n_[2]; @@ -37,7 +37,7 @@ poisson::poisson( std::shared_ptr > pptr ) klo = (int) 1.0*(ipslow)*nslow/npslow; khi = (int) 1.0*(ipslow+1)*nslow/npslow - 1; - //define local coords for wavenumbers + //Define local coords for wavenumber calculations inxlo = static_cast (1.0 * ipfast * nfast / npfast); inxhi = static_cast (1.0 * (ipfast+1) * nfast / npfast) - 1; @@ -49,6 +49,7 @@ poisson::poisson( std::shared_ptr > pptr ) nfft_in = (inxhi-inxlo+1) * (inyhi-inylo+1) * (inzhi-inzlo+1); +//Indexing for MPI running { for (int iproc = 0; iproc < parallelCom::worldSize(); iproc++) { if (parallelCom::myProcNo() != iproc) continue; @@ -78,21 +79,71 @@ poisson::poisson( std::shared_ptr > pptr ) if (parallelCom::myProcNo() < parallelCom::worldSize()) MPI_Send(&tmp,0,MPI_INT,parallelCom::myProcNo()+1,0,world); } } - //sanity check for fft coordinates - printf("Assigned local start coordinates (%d,%d,%d) on proc %d \n", - localStart_[0],localStart_[1],localStart_[2],parallelCom::myProcNo()); - printf("local range ijk (%d,%d,%d,%d,%d,%d)\n", inxlo,inxhi,inylo,inyhi,inzlo,inzhi); + //Sanity check for fft coordinates. Uncomment for debugging + //printf("Assigned local start coordinates (%d,%d,%d) on proc %d \n", + // localStart_[0],localStart_[1],localStart_[2],parallelCom::myProcNo()); + //printf("local range ijk (%d,%d,%d,%d,%d,%d)\n", inxlo,inxhi,inylo,inyhi,inzlo,inzhi); //initialize local fft domain + up_N = std::make_unique>(std::array ({inxlo, inylo, inzlo}) ,std::array ({inxhi, inyhi, inzhi})); N = up_N.get(); - //allocate memory for fft + up_NDCT = std::make_unique>(std::array ({inxlo, inylo, inzlo}) ,std::array ({inxhi, inyhi, inzhi})); + NDCT = up_NDCT.get(); + + //Construct and allocate memory for fft and dct + //HeFFTe::fftw_cos1 selects FFTW DCT TYPE REDFT00 + dct = std::make_unique>(*NDCT, *NDCT, world); + fft2d = std::make_unique>(*NDCT, *NDCT, world); + dct->nullify_executor(0); + dct->nullify_executor(2); + fft2d->nullify_executor(1); + + + fftx = std::make_unique>(*N, *N, world); + fftx->nullify_executor(1); + fftx->nullify_executor(2); + fft = std::make_unique>(*N, *N, world); + //plan 3d DCT-FFT-FFT + //perform 1d r2r_DCT + //translate domain from std::double to std::complex + //perform 2d r2c_DFT + //convolute and solve + //perform 2d c2r_IDFT + //translate domain from std::complex to std::double + //perform 1d r2r_DCT + //std::double cosphi + //std::complex dctphi + //std::complex dctphihat + //heffte::fft1d> + //heffte::fft2d> + + + //Create and allocate memory for FFT/DCT inputs & outputs. up_k2_ = std::make_unique((fft->size_inbox())); k2_ = up_k2_.get(); + + up_k2c_ = std::make_unique((2*dct->size_inbox())); + k2c_ = up_k2c_.get(); + + up_k2cy_ = std::make_unique((2*dct->size_inbox())); + k2cy_ = up_k2cy_.get(); + + up_dctworkspace_ = std::make_unique>(2*dct->size_workspace()); + dctworkspace_ = up_dctworkspace_.get(); + + up_cosphihat_ = std::make_unique>((2*dct->size_outbox())); + cosphihat_ = up_cosphihat_.get(); + + up_cosphi_ = std::make_unique>((2*dct->size_inbox())); + cosphi_ = up_cosphi_.get(); + + up_cosphicomplex_ = std::make_unique>>((2*fft2d->size_outbox())); + cosphicomplex_ = up_cosphicomplex_.get(); up_phi_ = std::make_unique>((fft->size_inbox())); phi_ = up_phi_.get(); @@ -100,7 +151,16 @@ poisson::poisson( std::shared_ptr > pptr ) up_phihat_ = std::make_unique>>((fft->size_outbox())); phihat_ = up_phihat_.get(); - wavenumbers_(); + up_psi_ = std::make_unique>((fftx->size_inbox())); + psi_ = up_psi_.get(); + + up_psihat_ = std::make_unique>>((fftx->size_outbox())); + psihat_ = up_psihat_.get(); + + + + wavenumbersFFT_(); + wavenumbersDCT_(); } @@ -120,11 +180,53 @@ void poisson::rhs( std::shared_ptr > f ) } } +void poisson::rhsDCT( std::shared_ptr > f ) +{ + + int l=0; + for( int k=(settings::m()/2); k<(pptr_->nk()-1-settings::m()/2); k++ ) + { + for( int j=(settings::m()/2); j<(pptr_->nj()-1-settings::m()/2); j++ ) + { + for( int i=(settings::m()/2); i<(pptr_->ni()-1-settings::m()/2); i++ ) + { + cosphi_[0][l] = f->operator()(i, j, k) ; + l++; + } + } + } + +} + +void poisson::x_FFT( std::shared_ptr > f ) +{ + int l=0; + for( int k=(settings::m()/2)+1; k<(settings::m()/2)+2; k++ ) + { + for( int j=(settings::m()/2); j<(pptr_->nj()-1-settings::m()/2); j++ ) + { + for( int i=(settings::m()/2); i<(pptr_->ni()-1-settings::m()/2); i++ ) + { + psi_[0][l] = f->operator()(i, j, k) ; + l++; + } + } + } + fftx->forward(psi_->data(), psihat_->data()); + + std::ofstream data("fftxdat.dat"); + for(size_t ll=0; ll<((fftx->size_outbox())/nk_); ll++) { + data<forward(phi_->data(), phihat_->data()); + //dealiasing +// scalar kk = (n_[0]*n_[1]*n_[2]) * 2.0/3.0; +// int l=0; for( int k=0; k(0.0, 0.0); } - else - { - phihat_[0][l] = std::complex(( -std::real(phihat_[0][l]) / k2_[l]) , -std::imag(phihat_[0][l]) / k2_[l]); - } + else + { + phihat_[0][l] = std::complex(( -std::real(phihat_[0][l]) / k2_[l]) , -std::imag(phihat_[0][l]) / k2_[l]); + } + l++; } - } + } } - fft->backward(phihat_->data(), phi_->data()); l=0; @@ -162,7 +264,77 @@ void poisson::solve() } } -void poisson::wavenumbers_() + +void poisson::solveDCT() +{ + dct->forward(cosphi_->data(), cosphihat_->data()); + + int l=0; + for( int k=0; kforward(cosphihat_->data(), cosphicomplex_->data()); + + l=0; + for( int k=0; k(0.0 , 0.0); + } + else + { + cosphicomplex_[0][l] =( std::complex(( -std::real(cosphicomplex_[0][l]) / k2c_[l]) , -std::imag(cosphicomplex_[0][l]) / k2c_[l] ) ); + } + l++; + } + } + } + + fft2d->backward(cosphicomplex_->data(), cosphihat_->data()); + + dct->backward(cosphihat_->data(), cosphi_->data(), dctworkspace_->data()); + + scalar pmax=0.0; + l=0; + for( int k=(settings::m()/2); k<(pptr_->nk()-1-settings::m()/2); k++ ) + { + for( int j=(settings::m()/2); j<(pptr_->nj()-1-settings::m()/2); j++ ) + { + for( int i=(settings::m()/2); i<(pptr_->ni()-1-settings::m()/2); i++ ) + { +// pmax = std::max( pmax, (cosphi_[0][l] / (2.0*n_[0]*n_[1]*(n_[2]-1) ))); + pptr_->operator()(i,j,k) = cosphi_[0][l] / (2*(n_[0])*(n_[1]-1)*n_[2]); + l++; + + } + } + } +// std::cout<<" "<size_inbox()<size_outbox()<size_outbox()<size_outbox()<size_outbox()<size_inbox()<mesh(); auto pi = tools::pi; @@ -223,10 +395,78 @@ void poisson::wavenumbers_() } } +void poisson::wavenumbersDCT_() +{ + Mesh& m = pptr_->mesh(); + auto pi = tools::pi; + + + scalar kx, ky, kz; //TODO: will be complex for non-central schemes + int l=0; + kx=0; + ky=0; + kz=0; + for( int k=0; k tools::eps ) + { + for( int ll=0; ll<(settings::m()/2); ll++ ) + { + kpx+=kx*settings::coef()[ll]/(ll+1)*std::sin((ll+1)*kx*m.dx())/(kx*m.dx()); + } + } + if(ky > tools::eps ) + { + for( int ll=0; ll<(settings::m()/2); ll++ ) + { + kpy+=ky*settings::coef()[ll]/(ll+1)*std::sin((ll+1)*ky*m.dy())/(ky*m.dy()); + } + } + if(kz > tools::eps ) + { + for( int ll=0; ll<(settings::m()/2); ll++ ) + { + kpz+=kz*settings::coef()[ll]/(ll+1)*std::sin((ll+1)*kz*m.dz())/(kz*m.dz()); + } + } + + k2cy_[l] = (kpy/(n_[1]-1)) * (kpy/(n_[1]-1) ); + k2c_[l] = kpx*kpx + kpz*kpz; + l++; + } + } + } +} + +//Destructor no longer needed, can remove poisson::~poisson() { free(N); fftw_free(phi_); + fftw_free(cosphi_); fftw_free(phihat_); } diff --git a/src/poisson/poisson.h b/src/poisson/poisson.h index 72dca51..78ed9c5 100644 --- a/src/poisson/poisson.h +++ b/src/poisson/poisson.h @@ -22,6 +22,7 @@ #include #include "heffte.h" +#include #include "Mesh/Mesh.h" #include "Field/Field.h" @@ -36,6 +37,9 @@ class poisson { protected: std::unique_ptr> fft; + std::unique_ptr> fft2d; + std::unique_ptr> dct; + std::unique_ptr> fftx; std::shared_ptr > pptr_; int fftsize; unsigned ni_, nj_, nk_; @@ -61,22 +65,53 @@ class poisson std::unique_ptr>> up_phihat_; std::vector>* phihat_; + std::unique_ptr> up_cosphi_; + std::vector* cosphi_; + + std::unique_ptr>> up_cosphicomplex_; + std::vector>* cosphicomplex_; + + std::unique_ptr> up_cosphihat_; + std::vector* cosphihat_; + + std::unique_ptr> up_dctworkspace_; + std::vector* dctworkspace_; + std::unique_ptr> up_phi_; std::vector* phi_; + std::unique_ptr>> up_psihat_; + std::vector>* psihat_; + + std::unique_ptr> up_psi_; + std::vector* psi_; + + std::unique_ptr up_k2cy_; + double* k2cy_; + std::unique_ptr up_k2_; double* k2_; - + + std::unique_ptr up_k2c_; + double* k2c_; + heffte::box3d<>* N; std::unique_ptr> up_N; + heffte::box3d<>* NDCT; + std::unique_ptr> up_NDCT; + public: + poisson( std::shared_ptr > ); ~poisson(); - void wavenumbers_(); - + void wavenumbersFFT_(); + void wavenumbersDCT_(); + void x_FFT( std::shared_ptr > f ); void rhs(std::shared_ptr > f ); + void rhsDCT(std::shared_ptr > f ); void solve(); + void solveDCT(); }; #endif diff --git a/src/reuseTmp/reuseTmp.hxx b/src/reuseTmp/reuseTmp.hxx index be178e9..4b8b1aa 100644 --- a/src/reuseTmp/reuseTmp.hxx +++ b/src/reuseTmp/reuseTmp.hxx @@ -32,7 +32,8 @@ reuseTmp::reuseTmp for( int i=0; i<=20; i++ ) { assert( i!=20 ); - if( tmp_[i].unique() || tmp_[i] == NULL) +// if( tmp_[i].unique() || tmp_[i] == NULL) + if( tmp_[i].use_count() == 1 || tmp_[i] == NULL) { if( tmp_[i] == NULL ) { From 92a1182d4d070c5467c68ff4d8ec7fcf12b6a294 Mon Sep 17 00:00:00 2001 From: u43070nw Date: Tue, 28 Jan 2025 17:26:42 +0000 Subject: [PATCH 3/3] Added HDF5 file writing in parallel. --- .DS_Store | Bin 6148 -> 0 bytes app/.DS_Store | Bin 6148 -> 0 bytes app/mhdjet/mhdjet.cpp | 4 +--- 3 files changed, 1 insertion(+), 3 deletions(-) delete mode 100644 .DS_Store delete mode 100644 app/.DS_Store diff --git a/.DS_Store b/.DS_Store deleted file mode 100644 index fd28e08626419453760300e8908c97e8b9d4f55a..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 6148 zcmeHKOH0E*5Z-O8-BN@e6nb3nTCfivh?fxSA26Z^m736^!8BW%)Ci@Jv;HA}iND90 z-3?gUgC`L?1DkJl9=kgqWFL$%?#_lC#u|(<0S%F(vPRIn)>X2>h+K~`VG$d#B*>D< znt}eJ3D+)Q&9BU7WA<(RL8JokeFT#zOwxAylb6bst?jC)in_S>o>b#9A$&X?zxInKa59Hoe&O35OQ}D#zU34Y8DR@o$DEgsEOL3(VWkF-J_Q5 zopcs0IX~^TTk^1Xyjax4-u}Va<={CQ$Lh_H$$@t*I~EIg2W6$M7jGIzDt-iCkypeM z5(C5lF+dD#4g>ZK5cSO|nJOg)h=Ct6fct}lhUi$#4eG4}K6rgbe+3Z*bbL!7N{f!g z+#q;BxJd;xsoXv>xJd`Qv~iBb+@MKkT+a;S*qMvR3)i!QUFvYg9fR~G28e-W2Fj+} z!1MnCewnq8{N)tt5d*})KVyJb`d;6KqU_naQXZbQ650ba6pYJJ0ResK5&#|CN1EDc c{1SDDb1dctaTK)cbU?ZYXhP^i4EzEEU;6${9smFU diff --git a/app/.DS_Store b/app/.DS_Store deleted file mode 100644 index d3a5c486434b2765aeb4cf95c8447bc1844da532..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 6148 zcmeHKQA@)x5WdW*Eu!o}!N-8F12=UJ@uf_C7c2UpGFv;eSQ}|O_b>*1)<5Jg@%MO_ zWa60mB8ZGTxO|t(T|&N;Tmt~2(~p_}H2`o>2`er(UkHtpPD#aj2!);_g&YdVpbs}w z(d_t-4A9xF;EDC{#1wo!zXFaSfGiy3k71%hjC&0qlQ=ILjWB&dP@*Wi=M^OAda92lZt3kg?(ZOla79A<2;L{L6Z(buZ;89m4$tw2)#P`r49$- z8RV83U jLziNVrBb|(ss;U$3`EajX%Iaq{3D=g;D#CaQwH7vIT=x6 diff --git a/app/mhdjet/mhdjet.cpp b/app/mhdjet/mhdjet.cpp index ed42e26..b9949a2 100644 --- a/app/mhdjet/mhdjet.cpp +++ b/app/mhdjet/mhdjet.cpp @@ -183,7 +183,6 @@ int main(int argc, char* argv[]) { for( int k=settings::m()/2; k