Skip to content

Commit 67f88ef

Browse files
committed
Added utility to write plotflile to hdf5. Added script to plot hdf5 file using matplotlib. Added run script to automate converting all the plotfiles in a directory and plotting them.
1 parent d8db933 commit 67f88ef

File tree

5 files changed

+258
-0
lines changed

5 files changed

+258
-0
lines changed

swm_AMReX/plotting_utils/GNUmakefile

+34
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,34 @@
1+
# AMREX_HOME defines the directory in which we will find all the AMReX code.
2+
AMREX_HOME ?= PROVIDE_YOUR_OWN_PATH_TO_AMREX_HERE
3+
HDF5_HOME ?= PROVIDE_YOUR_OWN_PATH_TO_HDF5_HERE
4+
5+
DEBUG = FALSE
6+
USE_MPI = TRUE # Dont really nedd MPI for the program but I am using an hdf5 that needs MPI
7+
USE_OMP = FALSE
8+
COMP = gnu
9+
DIM = 2
10+
11+
USE_HDF5 = TRUE
12+
13+
include $(AMREX_HOME)/Tools/GNUMake/Make.defs
14+
15+
include ./Make.package
16+
17+
## Set the target name explicitly
18+
#TARGET = plotfile2hdf5
19+
## Set the target name explicitly
20+
#EXEC = plotfile2hdf5
21+
#
22+
#executable = plotfile2hdf5.exe
23+
24+
include $(AMREX_HOME)/Src/Base/Make.package
25+
26+
vpath %.c : . $(vpathdir)
27+
vpath %.h : . $(vpathdir)
28+
vpath %.cpp : . $(vpathdir)
29+
vpath %.H : . $(vpathdir)
30+
vpath %.F : . $(vpathdir)
31+
vpath %.f : . $(vpathdir)
32+
vpath %.f90 : . $(vpathdir)
33+
34+
include $(AMREX_HOME)/Tools/GNUMake/Make.rules

swm_AMReX/plotting_utils/Make.package

+5
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,5 @@
1+
#CEXE_sources += plotfile2hdf5.cpp
2+
CEXE_sources += main.cpp
3+
4+
# Executable name
5+
#EXE = plotfile2hdf5
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,63 @@
1+
#!/bin/bash
2+
3+
# This script is used to convert all AMReX plotfiles in a specified directory to HDF5 format
4+
# and then plot the hdf5 files using matplotlib.
5+
6+
###############################################################################
7+
# User Input
8+
###############################################################################
9+
plotfile2hdf5_exe="$SWM_AMREX_ROOT"/plotting_utils/main2d.gnu.MPI.ex
10+
11+
# Input directory argument must be provided
12+
if [ -z "$1" ]; then
13+
echo "Usage: $0 <directory with amrex plotfiles of the form plt_*>"
14+
exit 1
15+
fi
16+
17+
# Set the directory to the provided argument
18+
input_directory="$1"
19+
20+
# Remove trailing slash from input_directory if it exists... so the pattern in the for loop below works
21+
input_directory="${input_directory%/}"
22+
23+
###############################################################################
24+
25+
set -u
26+
set -e
27+
28+
# Convenience function to print 80 character wide banners with centered text
29+
print_banner() {
30+
local message="$1"
31+
local border_char="#"
32+
local width=80
33+
local padding=$(( (width - ${#message} - 2) / 2 ))
34+
local border=$(printf "%-${width}s" | tr ' ' "$border_char")
35+
36+
echo -e "\n$border"
37+
printf "%*s %s %*s\n" $padding "" "$message" $padding ""
38+
echo -e "$border\n"
39+
}
40+
41+
42+
# Loop over all files that start with plt_ followed by a series of numbers in the specified directory
43+
for infile in "$input_directory"/plt[0-9]*; do
44+
45+
# Skip files that contain .old or .h5
46+
if [[ "$infile" == *".old"* || "$infile" == *".h5"* ]]; then
47+
continue
48+
fi
49+
50+
print_banner "Converting $infile to HDF5 and Plotting"
51+
52+
# Plotfiles are actually directories. Check that infile is a directory
53+
if [[ -d "$infile" ]]; then
54+
# Set the outfile name to the same as infile but with .h5 extension
55+
hdf5_file="${infile}.h5"
56+
57+
# Run the program with the infile and outfile arguments
58+
"$plotfile2hdf5_exe" infile="$infile" outfile="$hdf5_file"
59+
fi
60+
61+
python "$SWM_AMREX_ROOT"/plotting_utils/plot_hdf5.py "$hdf5_file" --output_dir "$input_directory"
62+
63+
done

swm_AMReX/plotting_utils/main.cpp

+101
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,101 @@
1+
#include <string>
2+
3+
#include <AMReX.H>
4+
#include <AMReX_ParmParse.H>
5+
#include <AMReX_PlotFileUtil.H>
6+
#include <AMReX_MultiFab.H>
7+
#include <hdf5.h>
8+
9+
void WriteMultiFabToHDF5(const std::string& plotfile, const std::string& hdf5file) {
10+
11+
amrex::PlotFileData plotfile_data(plotfile);
12+
amrex::Vector<std::string> varnames = plotfile_data.varNames();
13+
double time = plotfile_data.time();
14+
int time_step = plotfile_data.levelStep(0);
15+
16+
amrex::MultiFab mf;
17+
18+
amrex::VisMF::Read(mf, plotfile+"/Level_0/Cell");
19+
20+
amrex::BoxArray ba = mf.boxArray();
21+
22+
// Expecting a cell centered box with low and high index bounds {0,0} to {nx-1,ny-1}
23+
amrex::Box minimal_box = ba.minimalBox();
24+
AMREX_ASSERT(minimal_box.smallEnd(0) == 0);
25+
AMREX_ASSERT(minimal_box.smallEnd(1) == 0);
26+
27+
hid_t file_id = H5Fcreate(hdf5file.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
28+
29+
{
30+
hid_t attr_space_id = H5Screate(H5S_SCALAR);
31+
hid_t attr_id = H5Acreate(file_id, "time", H5T_NATIVE_DOUBLE, attr_space_id, H5P_DEFAULT, H5P_DEFAULT);
32+
H5Awrite(attr_id, H5T_NATIVE_DOUBLE, &time);
33+
H5Aclose(attr_id);
34+
H5Sclose(attr_space_id);
35+
}
36+
37+
{
38+
hid_t attr_space_id = H5Screate(H5S_SCALAR);
39+
hid_t attr_id = H5Acreate(file_id, "time_step", H5T_NATIVE_INT, attr_space_id, H5P_DEFAULT, H5P_DEFAULT);
40+
H5Awrite(attr_id, H5T_NATIVE_INT, &time_step);
41+
H5Aclose(attr_id);
42+
H5Sclose(attr_space_id);
43+
}
44+
45+
// Get the dimensions of the MultiFab
46+
hsize_t dims[2] = {static_cast<hsize_t>(minimal_box.length(0)), static_cast<hsize_t>(minimal_box.length(1))};
47+
48+
// Iterate over the components of the MultiFab
49+
for (int component_idx = 0; component_idx < mf.nComp(); ++component_idx) {
50+
// Create a dataset for this component
51+
hid_t dataspace_id = H5Screate_simple(2, dims, NULL);
52+
hid_t dataset_id = H5Dcreate(file_id, varnames[component_idx].c_str(), H5T_NATIVE_DOUBLE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
53+
54+
// Write the data to the dataset
55+
for (amrex::MFIter mfi(mf); mfi.isValid(); ++mfi) {
56+
const amrex::Box& bx = mfi.validbox();
57+
const amrex::Array4<const amrex::Real> &fab = mf.array(mfi);
58+
std::vector<double> data(bx.numPts());
59+
int idx = 0;
60+
for (int j = bx.smallEnd(1); j <= bx.bigEnd(1); ++j) {
61+
for (int i = bx.smallEnd(0); i <= bx.bigEnd(0); ++i) {
62+
data[idx++] = fab(i, j, 0, component_idx);
63+
}
64+
}
65+
H5Dwrite(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, data.data());
66+
}
67+
68+
H5Dclose(dataset_id);
69+
H5Sclose(dataspace_id);
70+
}
71+
72+
H5Fclose(file_id);
73+
}
74+
75+
int main(int argc, char* argv[]) {
76+
amrex::Initialize(argc, argv);
77+
78+
// Read in parameters from inputs file
79+
amrex::ParmParse pp;
80+
81+
// Read in plotfile name
82+
std::string plotfile;
83+
pp.query("infile", plotfile);
84+
if (plotfile.empty()) {
85+
amrex::Abort("You must specify `infile'");
86+
}
87+
88+
// Read in optional output hdf5 file name
89+
// Output hdf5 file (default to input plotfile name with .h5 extension)
90+
std::string hdf5file = plotfile + ".h5";
91+
pp.query("outfile", hdf5file);
92+
93+
amrex::Print() << "Input Plotfile: " << plotfile << std::endl;
94+
95+
WriteMultiFabToHDF5(plotfile, hdf5file);
96+
97+
amrex::Print() << "Output HDF5 file: " << hdf5file << std::endl;
98+
99+
amrex::Finalize();
100+
return 0;
101+
}

swm_AMReX/plotting_utils/plot_hdf5.py

+55
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,55 @@
1+
import argparse
2+
import os
3+
import h5py
4+
import numpy as np
5+
import matplotlib.pyplot as plt
6+
7+
if __name__ == "__main__":
8+
9+
parser = argparse.ArgumentParser(description='Read datasets from an HDF5 file. Save a plot of each dataset as a PNG file.')
10+
parser.add_argument('hdf5file', type=str, help='Path to the HDF5 file to read.')
11+
parser.add_argument('--output_dir', type=str, default=os.getcwd(), help='Directory to save the output PNG files to.')
12+
13+
args = parser.parse_args()
14+
15+
title_fontsize = 18
16+
label_fontsize = 12
17+
18+
with h5py.File(args.hdf5file, 'r') as f:
19+
# Read the time attribute
20+
time_step = f.attrs['time_step']
21+
time = f.attrs['time']
22+
23+
# Iterate over all datasets in the file
24+
for dataset_name in f:
25+
dataset = f[dataset_name]
26+
data = dataset[:]
27+
data_2d = np.array(data).reshape(dataset.shape)
28+
29+
plt.figure()
30+
plt.imshow(data_2d, aspect='auto', cmap='viridis')
31+
plt.colorbar()
32+
plt.title(f"${dataset_name}$", fontsize=title_fontsize)
33+
plt.xlabel('$x$', fontsize=label_fontsize)
34+
plt.ylabel('$y$', fontsize=label_fontsize)
35+
36+
# Print time step in the top left corner
37+
plt.gcf().text(0, 1.01,
38+
f"Time Step: {time_step}",
39+
transform=plt.gca().transAxes,
40+
fontsize=label_fontsize,
41+
verticalalignment='bottom',
42+
horizontalalignment='left')
43+
44+
# Print time in the top right corner
45+
plt.gcf().text(1.00, 1.01,
46+
f"Time: {time}",
47+
transform=plt.gca().transAxes,
48+
fontsize=label_fontsize,
49+
verticalalignment='bottom',
50+
horizontalalignment='right')
51+
52+
output_file = os.path.join(args.output_dir, f"{dataset_name}_{time_step:05d}.png")
53+
plt.savefig(output_file)
54+
55+
print(f"Saved plot to {output_file}")

0 commit comments

Comments
 (0)