Skip to content

Added basic HDF5 support #3

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
ndw108 wants to merge 3 commits into
base: main
Choose a base branch
from
Open

Added basic HDF5 support #3

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
58bc3f1
Change to HeFFTe & FFTW3 backend
Mar 1, 2024
3c6b179
Added HDF5 file writing in parallel.
Jan 28, 2025
92a1182
Added HDF5 file writing in parallel.
Jan 28, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
20 changes: 13 additions & 7 deletions CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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...
Expand All @@ -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)
Expand All @@ -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)
Expand All @@ -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)
Expand Down
17 changes: 13 additions & 4 deletions app/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,15 +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} ${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)

12 changes: 6 additions & 6 deletions app/OT/BEqn.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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() )
{
Expand Down
54 changes: 43 additions & 11 deletions app/OT_Poisson/OT.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.001, 2.01, 10 ); //args: dt, endT, write interval / steps

const scalar pi = tools::pi;
parallelCom::decompose( settings::zoneName()+"/"+"mesh" );
Expand All @@ -45,26 +45,51 @@ int main(int argc, char* argv[])
std::shared_ptr<Field<scalar> > pB_ptr( std::make_shared<Field<scalar> >( 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);

std::default_random_engine eng( parallelCom::myProcNo() );
std::uniform_real_distribution<double> dist(-1.0,1.0);

/*
//initial conditions
for( int i=settings::m()/2; i<mesh.ni()-settings::m()/2; i++ )
{
for( int j=settings::m()/2; j<mesh.nj()-settings::m()/2; j++ )
{
for( int k=settings::m()/2; k<mesh.nk()-settings::m()/2; k++ )
{
scalar x = -pi+(i-settings::m()/2)*mesh.dx()+mesh.origin().x();
scalar y = -pi+(j-settings::m()/2)*mesh.dy()+mesh.origin().y();
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);
B(i, j, k) = vector(0.2, 0, 0);

}
}
}
*/

//initial conditions
for( int i=settings::m()/2; i<mesh.ni()-settings::m()/2; i++ )
{
for( int j=settings::m()/2; j<mesh.nj()-settings::m()/2; j++ )
{
for( int k=settings::m()/2; k<mesh.nk()-settings::m()/2; k++ )
{
//OT Vortex initial conditions
scalar x = (i-settings::m()/2)*mesh.dx()+mesh.origin().x();
scalar y = (j-settings::m()/2)*mesh.dy()+mesh.origin().y();
scalar z = (k-settings::m()/2)*mesh.dz()+mesh.origin().z();

U(i, j, k) = vector((-2.0)*std::sin(y),2.0*std::sin(x),0.0);
B(i, j, k) = 0.818*vector( ((-2.0)*std::sin(y*2.0))+std::sin(z),((2.0*std::sin(x))+std::sin(z)),( std::sin(x)+std::sin(y) ));
}
B(i, j, k) = 0.8*vector( ((-2.0)*std::sin(2.0*y))+std::sin(z),((2.0*std::sin(x))+std::sin(z)),( std::sin(x)+std::sin(y) ));

}
}
}

Expand All @@ -75,6 +100,7 @@ int main(int argc, char* argv[])

std::ofstream data( settings::zoneName()+"/"+"dat.dat");
scalar EkOld = 0.0;
scalar EmOld = 0.0;

boost::timer::cpu_timer timer;

Expand All @@ -84,7 +110,7 @@ int main(int argc, char* argv[])

time++;

#include "RKLoop.H"
#include "UEqn.H"
#include "BEqn.H"

mesh.write(settings::zoneName()+"/data");
Expand All @@ -94,20 +120,21 @@ int main(int argc, char* argv[])
J.write(settings::zoneName()+"/data", "J");
pB.write(settings::zoneName()+"/data", "pB");

if( time.writelog() )
{
// if( time.writelog() )
// {
if( parallelCom::master() )
{
std::cout << "Step: " << time.timeStep() << ". Time: " << time.curTime() << std::endl;
std::cout << "Elapsed wall time for timestep: " << timer.elapsed().wall / 1e9 - start <<std::endl;
}

tools::CFL( U, mesh );
}
// }

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;
Expand All @@ -122,6 +149,7 @@ int main(int argc, char* argv[])
Ek += 0.5 * (U(i, j, k).x() * U(i, j, k).x() + U(i, j, k).y() * U(i, j, k).y() + U(i, j, k).z() * U(i, j, k).z() );
Em += 0.5 * (B(i, j, k).x() * B(i, j, k).x() + B(i, j, k).y() * B(i, j, k).y() + B(i, j, k).z() * B(i, j, k).z() );
Jmax = std::max( Jmax, sqrt(J(i, j, k).x() * J(i, j, k).x() + J(i, j, k).y() * J(i, j, k).y() + J(i, j, k).z() * J(i, j, k).z()));
Umax = std::max( Umax, sqrt(U(i, j, k).x() * U(i, j, k).x() + U(i, j, k).y() * U(i, j, k).y() + U(i, j, k).z() * U(i, j, k).z()));
omega += sqrt(omegav(i, j, k).x() * omegav(i, j, k).x() + omegav(i, j, k).y() * omegav(i, j, k).y() + omegav(i, j, k).z() * omegav(i, j, k).z());
Jav += sqrt(J(i, j, k).x() * J(i, j, k).x() + J(i, j, k).y() * J(i, j, k).y() + J(i, j, k).z() * J(i, j, k).z());
n++;
Expand All @@ -131,6 +159,7 @@ int main(int argc, char* argv[])

reduce( Ek, plusOp<scalar>() );
reduce( Jmax, maxOp<scalar>() );
reduce( Umax, maxOp<scalar>() );
reduce( Em, plusOp<scalar>() );
reduce( Jav, plusOp<scalar>() );
reduce( omega, plusOp<scalar>() );
Expand All @@ -141,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<<time.curTime()<<" "<<std::setprecision(15)<<Ek<<" "<<std::setprecision(15)<<Em<<" "<<std::setprecision(15)<<Jav<<" "<<std::setprecision(15)<<Jmax<<" "<<std::setprecision(15)<<omega<<" "<<std::setprecision(15)<<epsilon<<" "<<std::setprecision(15)<<-(Ek-EkOld)/time.dt()<<std::endl;
data<<time.curTime()<<" "<<std::setprecision(15)<<Ek<<" "<<std::setprecision(15)<<Em<<" "<<std::setprecision(15)<<Jav<<" "<<std::setprecision(15)<<Jmax<<" "<<std::setprecision(15)<<omega<<" "<<std::setprecision(15)<<epsilon<<" "<<std::setprecision(15)<<-(Ek-EkOld)/time.dt()<<std::setprecision(15)<<" "<<(Em-EmOld)/time.dt()<<std::setprecision(15)<<" "<<Umax<<std::endl;
}

EkOld = Ek;
EmOld = Em;

}

Expand Down
7 changes: 4 additions & 3 deletions app/OT_Poisson/RKLoop.H → app/OT_Poisson/UEqn.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,14 +28,15 @@
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() )
Expand Down
2 changes: 1 addition & 1 deletion app/TG/TG.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@
int main(int argc, char* argv[])
{
settings::process( argc, argv );
Time time( 0.001, 20, 1000 ); //args: dt, endT, write interval / steps
Time time( 0.01, 6.01, 100 ); //args: dt, endT, write interval / steps

const scalar pi = 3.1415926536;
parallelCom::decompose( settings::zoneName()+"/"+"mesh" );
Expand Down
13 changes: 8 additions & 5 deletions app/TG_Piosson/RKLoop.H
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
{

std::array<std::shared_ptr<Field<vector> >, 4> dU;
for( int i=0; i<4; i++ )
{
Expand All @@ -25,16 +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() )
{
Expand Down
9 changes: 5 additions & 4 deletions app/TG_Piosson/TG.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@
int main(int argc, char* argv[])
{
settings::process( argc, argv );
Time time( 0.01, 20, 1.0 ); //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" );
Expand Down Expand Up @@ -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);
}
}
}
Expand All @@ -79,16 +80,16 @@ int main(int argc, char* argv[])
p.write(settings::zoneName()+"/data", "p");


if( time.writelog() )
{
// if( time.writelog() )
// {
if( parallelCom::master() )
{
std::cout << "Step: " << time.timeStep() << ". Time: " << time.curTime() << std::endl;
std::cout << "Elapsed wall time for timestep: " << timer.elapsed().wall / 1e9 - start <<std::endl;
}

tools::CFL( U, mesh );
}
// }

scalar Ek=0.0;
int n=0;
Expand Down
56 changes: 56 additions & 0 deletions app/mhdjet/BEqn.H
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,56 @@
{
std::array<std::shared_ptr<Field<vector> >, 4> dB;
for( int i=0; i<4; i++ )
{
reuseTmp<vector> dBt( mesh );
dB[i]= dBt();
}


reuseTmp<vector> Boldt( mesh );
std::shared_ptr<Field<vector> > 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<<std::endl;
}
}

}
Loading