Skip to content

Commit

Permalink
Merge pull request #195 from geodynamics/developer
Browse files Browse the repository at this point in the history 
Merge developer to master
  • Loading branch information
@andersp
andersp authored Jul 19, 2023
2 parents 77a0d09 + 260d714 commit 13e6d43
Show file tree
Hide file tree
Showing 10 changed files with 276 additions and 167 deletions.
42 changes: 42 additions & 0 deletions configs/Dockerfile
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,42 @@
FROM ubuntu:20.04

WORKDIR /home

# build environment
RUN apt-get update && \
DEBIAN_FRONTEND='noninteractive' \
DEBCONF_NONINTERACTIVE_SEEN='true' \
apt-get install --yes wget build-essential git cmake gfortran openmpi-bin && \
\
\
# install hdf5 1.12.1, and remain packages \
wget https://support.hdfgroup.org/ftp/HDF5/releases/hdf5-1.12/hdf5-1.12.1/src/hdf5-1.12.1.tar.gz && \
tar xzf hdf5-1.12.1.tar.gz && rm hdf5-1.12.1.tar.gz && cd hdf5-1.12.1 && ./configure --prefix=/usr/local/hdf5 --enable-parallel --enable-build-mode=production && \
make -j4 && make install && cd .. && \
apt-get install --yes libproj-dev libhdf5-mpi-dev libfftw3-dev libfftw3-mpi-dev \
liblapack-dev python3-h5py libcurl3-dev && \
\
\
# install zfp and H5Z-ZFP \
git -c http.sslVerify=false clone https://github.com/LLNL/zfp.git && cd zfp/ && \
sed -i 's/# DEFS += -DBIT_STREAM_WORD_TYPE=uint8/DEFS += -DBIT_STREAM_WORD_TYPE=uint8/g' Config && \
make && cd .. && \
git -c http.sslVerify=false clone https://github.com/LLNL/H5Z-ZFP.git && cd H5Z-ZFP/ && \
make FC= CC=mpicc ZFP_HOME=/home/zfp HDF5_HOME=/usr/local/hdf5 all install && \
cd ..

# changing BUILD_TIME will break the cache here so that sw4 can be built each time
ARG BUILD_TIME=unknown

# build sw4
RUN git -c http.sslVerify=false clone https://github.com/geodynamics/sw4.git && \
cd sw4/ && git switch developer && rm -r .git/ && \
make sw4 CXX=mpicxx FC=gfortran debug=no proj=yes hdf5=yes fftw=yes zfp=yes prec=double EXTRA_LINK_FLAGS="-L/usr/lib64 -lgfortran -lhdf5 -llapack" HDF5ROOT=/usr/local/hdf5 FFTWHOME=/usr/lib/x86_64-linux-gnu ZFPROOT=/home/zfp H5ZROOT=/home/H5Z-ZFP/install -j4

ENV PATH="$PATH:/home/sw4/optimize_mp" \
LD_LIBRARY_PATH="$LD_LIBRARY_PATH:/usr/local/hdf5/lib"

# run sw4
# I typically use:
# docker run --cap-add=SYS_PTRACE -it -v $(pwd):/home/model/ sw4:zfp
# to run the docker container in avoid of the CMA issue (https://github.com/open-mpi/ompi/issues/4948)
33 changes: 19 additions & 14 deletions src/ESSI3D.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -313,6 +313,7 @@ void ESSI3D::open_vel_file(int a_cycle, std::string& a_path, float_sw4 a_time,
Sarray& a_Z) {
bool debug = false;
/* debug = true; */
MPI_Comm comm = mEW->m_cartesian_communicator;

bool is_root = true;
int g = mEW->mNumberOfGrids - 1;
Expand Down Expand Up @@ -352,21 +353,13 @@ void ESSI3D::open_vel_file(int a_cycle, std::string& a_path, float_sw4 a_time,
if (!m_isRestart) {
m_hdf5helper->write_header(h, lonlat_origin, az, origin, a_cycle, a_time,
dt);

// Write z coodinates if necesito
if (mEW->topographyExists()) {
compute_image(a_Z, 0, 0);
if (m_precision == 4)
m_hdf5helper->write_topo(m_floatField[0]);
else if (m_precision == 8)
m_hdf5helper->write_topo(m_doubleField[0]);
}
}

if (debug)
cout << "Creating hdf5 velocity fields..." << endl;
}

MPI_Barrier(comm);

if (m_dumpInterval > 0) {
int nstep = (int)ceil(m_ntimestep / m_dumpInterval);
if (m_compressionMode > 0)
Expand All @@ -384,16 +377,27 @@ void ESSI3D::open_vel_file(int a_cycle, std::string& a_path, float_sw4 a_time,
m_bufferInterval);
}

MPI_Comm comm = mEW->m_cartesian_communicator;
MPI_Barrier(comm);

is_root = false;
m_hdf5helper->create_file(true, is_root);

m_hdf5_time += (MPI_Wtime() - hdf5_time);

m_nbufstep = 0;
m_fileOpen = true;

if (!m_isRestart) {
// Write z coodinates if has topo
if (mEW->topographyExists()) {
compute_image(a_Z, 0, 0);
if (m_precision == 4)
m_hdf5helper->write_topo(m_floatField[0]);
else if (m_precision == 8)
m_hdf5helper->write_topo(m_doubleField[0]);
}
}

m_hdf5_time += (MPI_Wtime() - hdf5_time);

return;
}

Expand All @@ -419,9 +423,10 @@ void ESSI3D::write_image_hdf5(int cycle, std::string& path, float_sw4 t,
/* debug = true; */

for (int i = 0; i < 3; i++) {
int nstep = (int)floor(m_ntimestep / m_dumpInterval);
compute_image(a_U[g], i, cycle);
if (cycle > 0 &&
(m_nbufstep == m_bufferInterval - 1 || cycle == m_ntimestep)) {
(m_nbufstep == m_bufferInterval - 1 || cycle == nstep)) {
if (m_precision == 4)
m_hdf5helper->write_vel((void*)m_floatField[i], i, cycle,
m_nbufstep + 1);
Expand Down
2 changes: 2 additions & 0 deletions src/ESSI3DHDF5.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -314,6 +314,8 @@ void ESSI3DHDF5::write_topo(void* window_array) {
ierr = H5Sclose(dataspace_id);
ierr = H5Dclose(dataset_id);

H5Fflush(m_file_id, H5F_SCOPE_LOCAL);

if (debug && (myRank == 0))
cerr << "Done writing hdf5 z coordinate: " << m_filename << endl;
#endif
Expand Down
4 changes: 3 additions & 1 deletion src/EW.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7238,7 +7238,6 @@ void EW::extractTopographyFromGMG( std::string a_topoFileName )
start_time = MPI_Wtime();
#ifdef USE_HDF5
int verbose = mVerbose;
std::string rname ="EW::extractTopographyFromGMG";
Sarray gridElev;
herr_t ierr;
hid_t file_id, dataset_id, datatype_id, group_id, dataspace_id;
Expand Down Expand Up @@ -8495,6 +8494,9 @@ void EW::set_to_zero_at_receiver( vector<Sarray> & a_U,
#pragma omp parallel for
for( int s=0 ; s < time_series.size() ; s++ )
{
if (!time_series[s]->myPoint())
continue;

int g = time_series[s]->m_grid0;
int i0= time_series[s]->m_i0;
int j0= time_series[s]->m_j0;
Expand Down
49 changes: 26 additions & 23 deletions src/Image.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1347,25 +1347,28 @@ void Image::writeImagePlane_2(int cycle, std::string &path, float_sw4 t )
H5Sclose(dset_space);
H5Dclose(dset);

if( mGridinfo == 1 )
char grid_name[16];
if( mGridinfo >= 1 )
{
int g=mEW->mNumberOfGrids-1;
int globalSizes[3] = {mEW->m_global_nx[g], mEW->m_global_ny[g], mEW->m_global_nz[g]} ;
if(mLocationType == Image::X)
globalSizes[0] = 1;
if (mLocationType == Image::Y)
globalSizes[1] = 1;
if (mLocationType == Image::Z)
globalSizes[2] = 1;

dims = globalSizes[0]*globalSizes[1]*globalSizes[2];
dset_space = H5Screate_simple(1, &dims, NULL);
/* cout << "Rank " << mEW->getRank() << " creating grid array with length " << dims << endl; */
dset = H5Dcreate(h5_fid, "grid", dtype, dset_space, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
if( dset < 0 )
cout << "ERROR: Image::writeImagePlane_2 could not create HDF5 grid dataset" << endl;
H5Sclose(dset_space);
H5Dclose(dset);
for (int g = mEW->mNumberOfCartesianGrids; g < mEW->mNumberOfGrids; g++) {
int globalSizes[3] = {mEW->m_global_nx[g], mEW->m_global_ny[g], mEW->m_global_nz[g]} ;
if(mLocationType == Image::X)
globalSizes[0] = 1;
if (mLocationType == Image::Y)
globalSizes[1] = 1;
if (mLocationType == Image::Z)
globalSizes[2] = 1;

dims = globalSizes[0]*globalSizes[1]*globalSizes[2];
dset_space = H5Screate_simple(1, &dims, NULL);
/* cout << "Rank " << mEW->getRank() << " creating grid array with length " << dims << endl; */
sprintf(grid_name, "grid%d", g);
dset = H5Dcreate(h5_fid, grid_name, dtype, dset_space, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
if( dset < 0 )
cout << "ERROR: Image::writeImagePlane_2 could not create HDF5 grid dataset" << endl;
H5Sclose(dset_space);
H5Dclose(dset);
}
}

delete [] grid_size;
Expand Down Expand Up @@ -1481,7 +1484,7 @@ void Image::writeImagePlane_2(int cycle, std::string &path, float_sw4 t )
}
} // End if ihavearray

if( mGridinfo == 1 )
if( mGridinfo >= 1 )
add_grid_to_file_hdf5( s.str().c_str(), iwrite, 0);
/* if( mGridinfo == 2 ) */
/* add_grid_filenames_to_file( s.str().c_str() ); */
Expand All @@ -1495,6 +1498,7 @@ void Image::writeImagePlane_2(int cycle, std::string &path, float_sw4 t )
void Image::add_grid_to_file_hdf5( const char* fname, bool iwrite, size_t offset )
{
bool ihavearray = plane_in_proc(m_gridPtIndex[0]);
char grid_name[16];
if( ihavearray )
{
#ifdef USE_HDF5
Expand All @@ -1512,17 +1516,16 @@ void Image::add_grid_to_file_hdf5( const char* fname, bool iwrite, size_t offset
if( mEW->usingParallelFS() )
MPI_Barrier( m_mpiComm_writers );

sprintf(grid_name, "grid%d", g);
if( m_double )
{
char dblStr[]="double";
m_pio[g]->write_array_hdf5(fname, NULL, "grid", 1, m_gridimage->m_doubleField[g], offset, dblStr );
offset += (globalSizes[0]*globalSizes[1]*globalSizes[2]*sizeof(double));
m_pio[g]->write_array_hdf5(fname, NULL, grid_name, 1, m_gridimage->m_doubleField[g], 0, dblStr );
}
else
{
char fltStr[]="float";
m_pio[g]->write_array_hdf5(fname, NULL, "grid", 1, m_gridimage->m_floatField[g], offset, fltStr );
offset += (globalSizes[0]*globalSizes[1]*globalSizes[2]*sizeof(float));
m_pio[g]->write_array_hdf5(fname, NULL, grid_name, 1, m_gridimage->m_floatField[g], 0, fltStr );
}
}
#endif
Expand Down
6 changes: 5 additions & 1 deletion src/Qspline.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -193,7 +193,11 @@ Qspline::Qspline( int npts, float_sw4* fun, float_sw4 tmin, float_sw4 dt, int bc
delete[] ipiv;

m_npts = npts;
m_polcof = new float_sw4[6*(npts-1)];
if (npts > 1)
m_polcof = new float_sw4[6*(npts-1)];
else
m_polcof = new float_sw4[6];

for( int i= 0 ; i < npts-1 ; i++ )
{
m_polcof[ 6*i] = fun[i];
Expand Down
43 changes: 28 additions & 15 deletions src/SfileOutput.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -351,8 +351,22 @@ void SfileOutput::compute_image( vector<Sarray>& a_U, vector<Sarray>& a_Rho,
{
int stH = mSampleH;
int stV = mSampleV;
vector<int> hh(mEW->mNumberOfGrids, 1);
double my_z, up_z, down_z, up_v, down_v;

// Calculate horizontal factor of each grid for topo data access
for (int i = mEW->mNumberOfGrids-2; i >= 0; i--) {
// No h increase from Cartesian to Curvilinear
if (i == mEW->mNumberOfCartesianGrids-1)
hh[i] = hh[i+1];
else
hh[i] = hh[i+1] * 2;
}

// debug
/* int myRank; */
/* MPI_Comm_rank( mEW->m_1d_communicator, &myRank); */

/* if (mMode== RHO && m_parallel_io[0]->proc_zero()) { */
/* for( int g=mEW->mNumberOfCartesianGrids ; g < mEW->mNumberOfGrids ; g++ ) { */
/* cout << "g="<< g << ", a_Z.m_kb="<< a_Z[g].m_kb << ", a_Z.m_ke=" << a_Z[g].m_ke << ", gz=" << mWindow[g][5] << endl; */
Expand Down Expand Up @@ -390,15 +404,15 @@ void SfileOutput::compute_image( vector<Sarray>& a_U, vector<Sarray>& a_Rho,

if( mMode == RHO || mMode == QP || mMode == QS ) { // these modes just copy the values straight from the array
if( m_double ) {
#pragma omp parallel for
/* #pragma omp parallel for */
for( int k=mWindow[g][4] ; k <= mWindow[g][5] ; k+=stV )
for( int j=mWindow[g][2] ; j <= mWindow[g][3] ; j+=stH )
for( int i=mWindow[g][0] ; i <= mWindow[g][1] ; i+=stH ) {
size_t ind = (k-mWindow[g][4])/stV+nkw*(j-mWindow[g][2])/stH+njkw*(i-mWindow[g][0])/stH;
if (g < mEW->mNumberOfCartesianGrids)
m_doubleField[g][ind] = (double)((*data1)(1,i,j,k));
else {
double z_kl = -mEW->mTopo(i,j,1);
double z_kl = a_Z[gz](i,j,kl);
double z_ku = a_Z[gz](i,j,ku);
my_z = z_kl + (z_ku - z_kl)*(k-1)/(double)(ku-1);
int t = 1;
Expand All @@ -415,16 +429,15 @@ void SfileOutput::compute_image( vector<Sarray>& a_U, vector<Sarray>& a_Rho,
}
}
else {
#pragma omp parallel for
/* #pragma omp parallel for */
for( int k=mWindow[g][4] ; k <= mWindow[g][5] ; k+=stV )
for( int j=mWindow[g][2] ; j <= mWindow[g][3] ; j+=stH )
for( int i=mWindow[g][0] ; i <= mWindow[g][1] ; i+=stH ) {
size_t ind = (k-mWindow[g][4])/stV+nkw*(j-mWindow[g][2])/stH+njkw*(i-mWindow[g][0])/stH;
if (g < mEW->mNumberOfCartesianGrids)
m_floatField[g][ind] = (float)((*data1)(1,i,j,k));
else {
// debug
double z_kl = -mEW->mTopo(i,j,1);
double z_kl = a_Z[gz](i,j,kl);
double z_ku = a_Z[gz](i,j,ku);
my_z = z_kl + (z_ku - z_kl)*(k-1)/(double)(ku-1);
int t = 1;
Expand All @@ -443,15 +456,15 @@ void SfileOutput::compute_image( vector<Sarray>& a_U, vector<Sarray>& a_Rho,
}
else if( mMode == P ) {
if( m_double ) {
#pragma omp parallel for
/* #pragma omp parallel for */
for( int k=mWindow[g][4] ; k <= mWindow[g][5] ; k+=stV )
for( int j=mWindow[g][2] ; j <= mWindow[g][3] ; j+=stH )
for( int i=mWindow[g][0] ; i <= mWindow[g][1] ; i+=stH ) {
size_t ind = (k-mWindow[g][4])/stV+nkw*(j-mWindow[g][2])/stH+njkw*(i-mWindow[g][0])/stH;
if (g < mEW->mNumberOfCartesianGrids)
m_doubleField[g][ind] = (double)sqrt((2*((*data2)(1,i,j,k)) +((*data3)(1,i,j,k)))/((*data1)(1,i,j,k)));
else {
double z_kl = -mEW->mTopo(i,j,1);
double z_kl = a_Z[gz](i,j,kl);
double z_ku = a_Z[gz](i,j,ku);
my_z = z_kl + (z_ku - z_kl)*(k-1)/(double)(ku-1);
int t = 1;
Expand All @@ -477,15 +490,15 @@ void SfileOutput::compute_image( vector<Sarray>& a_U, vector<Sarray>& a_Rho,
}
}
else {
#pragma omp parallel for
/* #pragma omp parallel for */
for( int k=mWindow[g][4] ; k <= mWindow[g][5] ; k+=stV )
for( int j=mWindow[g][2] ; j <= mWindow[g][3] ; j+=stH )
for( int i=mWindow[g][0] ; i <= mWindow[g][1] ; i+=stH ) {
size_t ind = (k-mWindow[g][4])/stV+nkw*(j-mWindow[g][2])/stH+njkw*(i-mWindow[g][0])/stH;
if (g < mEW->mNumberOfCartesianGrids)
m_floatField[g][ind] = (float)sqrt((2.0*((*data2)(1,i,j,k)) +((*data3)(1,i,j,k)))/((*data1)(1,i,j,k)));
else {
double z_kl = -mEW->mTopo(i,j,1);
double z_kl = a_Z[gz](i,j,kl);
double z_ku = a_Z[gz](i,j,ku);
my_z = z_kl + (z_ku - z_kl)*(k-1)/(double)(ku-1);
int t = 1;
Expand Down Expand Up @@ -520,15 +533,15 @@ void SfileOutput::compute_image( vector<Sarray>& a_U, vector<Sarray>& a_Rho,
}
else if( mMode == S ) {
if( m_double ) {
#pragma omp parallel for
/* #pragma omp parallel for */
for( int k=mWindow[g][4] ; k <= mWindow[g][5] ; k+=stV )
for( int j=mWindow[g][2] ; j <= mWindow[g][3] ; j+=stH )
for( int i=mWindow[g][0] ; i <= mWindow[g][1] ; i+=stH ) {
size_t ind = (k-mWindow[g][4])/stV+nkw*(j-mWindow[g][2])/stH+njkw*(i-mWindow[g][0])/stH;
if (g < mEW->mNumberOfCartesianGrids)
m_doubleField[g][ind] = (double)sqrt(((*data2)(1,i,j,k))/((*data1)(1,i,j,k)));
else {
double z_kl = -mEW->mTopo(i,j,1);
double z_kl = a_Z[gz](i,j,kl);
double z_ku = a_Z[gz](i,j,ku);
my_z = z_kl + (z_ku - z_kl)*(k-1)/(double)(ku-1);
int t = 1;
Expand All @@ -553,15 +566,15 @@ void SfileOutput::compute_image( vector<Sarray>& a_U, vector<Sarray>& a_Rho,
}
}
else {
#pragma omp parallel for
/* #pragma omp parallel for */
for( int k=mWindow[g][4] ; k <= mWindow[g][5] ; k+=stV )
for( int j=mWindow[g][2] ; j <= mWindow[g][3] ; j+=stH )
for( int i=mWindow[g][0] ; i <= mWindow[g][1] ; i+=stH ) {
size_t ind = (k-mWindow[g][4])/stV+nkw*(j-mWindow[g][2])/stH+njkw*(i-mWindow[g][0])/stH;
if (g < mEW->mNumberOfCartesianGrids)
m_floatField[g][ind] = (float)sqrt(((*data2)(1,i,j,k))/((*data1)(1,i,j,k)));
else {
double z_kl = -mEW->mTopo(i,j,1);
double z_kl = a_Z[gz](i,j,kl);
double z_ku = a_Z[gz](i,j,ku);
my_z = z_kl + (z_ku - z_kl)*(k-1)/(double)(ku-1);
int t = 1;
Expand Down Expand Up @@ -627,7 +640,7 @@ void SfileOutput::compute_image( vector<Sarray>& a_U, vector<Sarray>& a_Rho,
}
}
}
}
}// End for g < mEW->mNumberOfGrids
}

//-----------------------------------------------------------------------
Expand Down Expand Up @@ -894,7 +907,7 @@ void SfileOutput::write_image(const char *fname, std::vector<Sarray>& a_Z )
int nj = (int)(mWindow[real_g][3]-mWindow[real_g][2])/stH+1;
int nk = mWindow[real_g][5];
float* zfp = new float[npts];
#pragma omp parallel for
/* #pragma omp parallel for */
for( int j=mWindow[real_g][2] ; j <= mWindow[real_g][3] ; j+=stH )
for( int i=mWindow[real_g][0] ; i <= mWindow[real_g][1] ; i+=stH ) {
size_t ind = (size_t)(j-mWindow[real_g][2])/stH+nj*(i-mWindow[real_g][0])/stH;
Expand Down
Loading

0 comments on commit 13e6d43

Please sign in to comment.