diff --git a/tools/gridmanip/.gitignore b/tools/gridmanip/.gitignore new file mode 100644 index 0000000000000000000000000000000000000000..d923660a9a349e84cfdc76d0712097a3614c2239 --- /dev/null +++ b/tools/gridmanip/.gitignore @@ -0,0 +1,8 @@ +gridmanip +*.o +*.a +*~ +*.swp +*.bak +*.h5 +*.xmf diff --git a/tools/gridmanip/BlockProcessor_MPI.h b/tools/gridmanip/BlockProcessor_MPI.h new file mode 100644 index 0000000000000000000000000000000000000000..89d55ca8f8b4b107272bb65900116d362604a94f --- /dev/null +++ b/tools/gridmanip/BlockProcessor_MPI.h @@ -0,0 +1,69 @@ +// File : BlockProcessor_MPI.h +// Created : Sat May 04 2019 01:25:51 PM (+0200) +// Author : Fabian Wermelinger +// Description: Simple MPI block-processor +// Copyright 2019 ETH Zurich. All Rights Reserved. +#ifndef BLOCKPROCESSOR_MPI_H_VTNAQSNJ +#define BLOCKPROCESSOR_MPI_H_VTNAQSNJ + +#include "Cubism/SynchronizerMPI.h" + +#include +#include + +using namespace cubism; + +template +inline void +process(TKernel rhs, TGrid &grid, const Real t = 0.0, const bool record = false) +{ + TKernel myrhs = rhs; + + SynchronizerMPI &Synch = grid.sync(myrhs); + + std::vector avail0, avail1; + + const int nthreads = omp_get_max_threads(); + TLab *labs = new TLab[nthreads]; + for (int i = 0; i < nthreads; ++i) + labs[i].prepare(grid, Synch); + + MPI_Barrier(grid.getCartComm()); + + avail0 = Synch.avail_inner(); + BlockInfo *ary0 = &avail0.front(); +#pragma omp parallel num_threads(nthreads) + { + int tid = omp_get_thread_num(); + TLab &mylab = labs[tid]; + +#pragma omp for schedule(dynamic, 1) + for (size_t i = 0; i < avail0.size(); i++) { + mylab.load(ary0[i], t); + rhs(mylab, ary0[i], *(typename TGrid::BlockType *)ary0[i].ptrBlock); + } + } + + avail1 = Synch.avail_halo(); + BlockInfo *ary1 = &avail1.front(); +#pragma omp parallel num_threads(nthreads) + { + int tid = omp_get_thread_num(); + TLab &mylab = labs[tid]; + +#pragma omp for schedule(dynamic, 1) + for (size_t i = 0; i < avail1.size(); i++) { + mylab.load(ary1[i], t); + rhs(mylab, ary1[i], *(typename TGrid::BlockType *)ary1[i].ptrBlock); + } + } + + if (labs != NULL) { + delete[] labs; + labs = NULL; + } + + MPI_Barrier(grid.getCartComm()); +} + +#endif /* BLOCKPROCESSOR_MPI_H_VTNAQSNJ */ diff --git a/tools/gridmanip/Example/Makefile b/tools/gridmanip/Example/Makefile new file mode 100644 index 0000000000000000000000000000000000000000..2eb96e3ca7d051f9d811d1cef059a386396a2b46 --- /dev/null +++ b/tools/gridmanip/Example/Makefile @@ -0,0 +1,17 @@ +# File : Makefile +# Created : Wed Nov 07 2018 06:53:55 PM (+0100) +# Author : Fabian Wermelinger +# Description: Run test cases. You need to compile the gridmanip tool and +# create a copy or sym link in the directory your run this +# makefile. +# Copyright 2018 ETH Zurich. All Rights Reserved. + +.PHONY: all clean + +all: + ./gridmanip -conf smoother.conf + ./gridmanip -conf restrict.conf + ./gridmanip -conf prolong.conf + +clean: + rm -f *.h5 *.xmf diff --git a/tools/gridmanip/Example/data/nonuniform/data_000000-5.h5 b/tools/gridmanip/Example/data/nonuniform/data_000000-5.h5 new file mode 100644 index 0000000000000000000000000000000000000000..c07cae8ff0d126964ee773ddb006c2234617f939 Binary files /dev/null and b/tools/gridmanip/Example/data/nonuniform/data_000000-5.h5 differ diff --git a/tools/gridmanip/Example/data/nonuniform/data_000000-5.xmf b/tools/gridmanip/Example/data/nonuniform/data_000000-5.xmf new file mode 100644 index 0000000000000000000000000000000000000000..68b1be2209c36f52eb63775f388966a323c41320 --- /dev/null +++ b/tools/gridmanip/Example/data/nonuniform/data_000000-5.xmf @@ -0,0 +1,29 @@ + + + + + + + + diff --git a/tools/gridmanip/Example/data/uniform/data_000000-5.h5 b/tools/gridmanip/Example/data/uniform/data_000000-5.h5 new file mode 100644 index 0000000000000000000000000000000000000000..a04a036d491de00d4ec1785bc47529046595ee1e Binary files /dev/null and b/tools/gridmanip/Example/data/uniform/data_000000-5.h5 differ diff --git a/tools/gridmanip/Example/data/uniform/data_000000-5.xmf b/tools/gridmanip/Example/data/uniform/data_000000-5.xmf new file mode 100644 index 0000000000000000000000000000000000000000..68b1be2209c36f52eb63775f388966a323c41320 --- /dev/null +++ b/tools/gridmanip/Example/data/uniform/data_000000-5.xmf @@ -0,0 +1,29 @@ + + + + + + + + diff --git a/tools/gridmanip/Example/gridmanip b/tools/gridmanip/Example/gridmanip new file mode 120000 index 0000000000000000000000000000000000000000..79d138e1a6520808f6a93807c2f6652d3d59a509 --- /dev/null +++ b/tools/gridmanip/Example/gridmanip @@ -0,0 +1 @@ +../gridmanip \ No newline at end of file diff --git a/tools/gridmanip/Example/prolong.conf b/tools/gridmanip/Example/prolong.conf new file mode 100644 index 0000000000000000000000000000000000000000..96bd1ccfa968575ca4b5c94b829215bfd172ae03 --- /dev/null +++ b/tools/gridmanip/Example/prolong.conf @@ -0,0 +1,28 @@ +# File : prolong.conf +# Created : Wed Nov 07 2018 05:41:21 PM (+0100) +# Author : Fabian Wermelinger +# Description: Refine +# Copyright 2018 ETH Zurich. All Rights Reserved. +extent 1.0 + +xpesize 1 +ypesize 1 +zpesize 1 + +in_bpdx 2 +in_bpdy 2 +in_bpdz 2 + +out_bpdx 4 +out_bpdy 4 +out_bpdz 4 + +read_format cubism_h5 +# in_file data/nonuniform/data_000000-5 # non-uniform grid data +in_file data/uniform/data_000000-5 # uniform grid data + +save_format cubism_h5 +out_file prolongation + +operator ProlongHarten +smooth_iter 0 diff --git a/tools/gridmanip/Example/restrict.conf b/tools/gridmanip/Example/restrict.conf new file mode 100644 index 0000000000000000000000000000000000000000..04d9546415165047fab1e98b09be32e95db81fe4 --- /dev/null +++ b/tools/gridmanip/Example/restrict.conf @@ -0,0 +1,27 @@ +# File : restrict.conf +# Created : Wed Nov 07 2018 05:41:08 PM (+0100) +# Author : Fabian Wermelinger +# Description: Coarsen +# Copyright 2018 ETH Zurich. All Rights Reserved. +extent 1.0 + +xpesize 1 +ypesize 1 +zpesize 1 + +in_bpdx 2 +in_bpdy 2 +in_bpdz 2 + +out_bpdx 1 +out_bpdy 1 +out_bpdz 1 + +read_format cubism_h5 +# in_file data/nonuniform/data_000000-5 # non-uniform grid data +in_file data/uniform/data_000000-5 # uniform grid data + +save_format cubism_h5 +out_file restriction + +operator RestrictBlockAverage diff --git a/tools/gridmanip/Example/smoother.conf b/tools/gridmanip/Example/smoother.conf new file mode 100644 index 0000000000000000000000000000000000000000..0bb04f01370dbe0e8f98026e4f11c115e41a8d0f --- /dev/null +++ b/tools/gridmanip/Example/smoother.conf @@ -0,0 +1,28 @@ +# File : smoother.conf +# Created : Wed Nov 07 2018 05:39:18 PM (+0100) +# Author : Fabian Wermelinger +# Description: Smooth grid +# Copyright 2018 ETH Zurich. All Rights Reserved. +extent 1.0 + +xpesize 1 +ypesize 1 +zpesize 1 + +in_bpdx 2 +in_bpdy 2 +in_bpdz 2 + +out_bpdx 2 +out_bpdy 2 +out_bpdz 2 + +read_format cubism_h5 +# in_file data/nonuniform/data_000000-5 # non-uniform grid data +in_file data/uniform/data_000000-5 # uniform grid data + +save_format cubism_h5 +out_file smoothed + +operator Smoother +smooth_iter 200 diff --git a/tools/gridmanip/GridOperator.h b/tools/gridmanip/GridOperator.h new file mode 100644 index 0000000000000000000000000000000000000000..57eab30762044150c7d5e587ca060265aa47483a --- /dev/null +++ b/tools/gridmanip/GridOperator.h @@ -0,0 +1,30 @@ +// File : GridOperator.h +// Created : Mon Jul 10 2017 12:26:52 PM (+0200) +// Author : Fabian Wermelinger +// Description: Operator base class +// Copyright 2017 ETH Zurich. All Rights Reserved. +#ifndef GRIDOPERATOR_H_LBMHW8RU +#define GRIDOPERATOR_H_LBMHW8RU + +#include "Cubism/ArgumentParser.h" +#include "Types.h" + +using namespace cubism; + +template +class GridOperator +{ +public: + GridOperator(ArgumentParser& p) : m_parser(p) {} + + virtual ~GridOperator() = default; + + virtual void + operator()(const TGridIn &, TGridOut &, const bool verbose = true) = 0; + +protected: + ArgumentParser& m_parser; +}; + +#endif /* GRIDOPERATOR_H_LBMHW8RU */ + diff --git a/tools/gridmanip/Makefile b/tools/gridmanip/Makefile new file mode 100644 index 0000000000000000000000000000000000000000..cfbdacfb4190beddd9702ec19e76d6305a5eef95 --- /dev/null +++ b/tools/gridmanip/Makefile @@ -0,0 +1,62 @@ +SHELL := /bin/bash + +CC = mpic++ + +bs ?= 32 +align ?= 16 +omp ?= 1 +ap ?= float +config ?= release + +# HDF5 is required to build +hdf = 1 +hdf-inc ?= ${HDF5_ROOT}/include +hdf-lib ?= ${HDF5_ROOT}/lib + +GIT_COMMIT := -D_GIT_HASH_=\"$(shell git rev-parse --verify HEAD)\" +CPPFLAGS += $(GIT_COMMIT) -g -std=c++11 -D_HARTEN5_ +CPPFLAGS += $(extra) + +# compiler warnings +CPPFLAGS += -Wall -Wextra -Wfloat-equal -Wundef -Wcast-align +CPPFLAGS += -Wwrite-strings -Wmissing-declarations -Wredundant-decls +CPPFLAGS += -Wshadow -Woverloaded-virtual -Wuninitialized +CPPFLAGS += -Wpedantic -Wno-unused-parameter + +ifeq "$(omp)" "1" + CPPFLAGS += -fopenmp + OPTFLAGS += -fopenmp +endif + +ifeq "$(ap)" "float" + CPPFLAGS += -D_FLOAT_PRECISION_ +endif + +ifeq "$(config)" "release" + OPTFLAGS += -O3 -DNDEBUG +endif + +ifeq "$(hdf)" "1" + CPPFLAGS += -I$(hdf-inc) -DCUBISM_USE_HDF + LIBS += -L$(hdf-lib) -lhdf5 +endif + +CPPFLAGS += -DCUBISM_ALIGNMENT=$(align) -D_BLOCKSIZE_=$(bs) -D_BLOCKSIZEX_=$(bs) -D_BLOCKSIZEY_=$(bs) -D_BLOCKSIZEZ_=$(bs) +CPPFLAGS += -I. +CPPFLAGS += -I../../../Cubism/include + +.DEFAULT_GOAL := gridmanip + +OBJECTS = main.o \ + ../../../Cubism/src/ArgumentParser.o + +all: gridmanip + +gridmanip: $(OBJECTS) + $(CC) $(OPTFLAGS) $(extra) $^ -o $@ $(LIBS) + +%.o: %.cpp + $(CC) $(OPTFLAGS) $(CPPFLAGS) -c $^ -o $@ + +clean: + rm -f *.o gridmanip *~ diff --git a/tools/gridmanip/Prolongation/MPI_GridTransfer.h b/tools/gridmanip/Prolongation/MPI_GridTransfer.h new file mode 100644 index 0000000000000000000000000000000000000000..fdb0e7589d34a501f36af09a37048150946ebb45 --- /dev/null +++ b/tools/gridmanip/Prolongation/MPI_GridTransfer.h @@ -0,0 +1,346 @@ +/* File: MPI_GridTransfer.h */ +/* Date: August 2015 */ +/* Author: Ursula Rasthofer */ +/* Tag: Transfer solution fields from coarse to fine grid */ +/* Copyright 2015 ETH Zurich. All Rights Reserved. */ +#ifndef MPI_GRIDTRANSFER_H_JCMHPZQU +#define MPI_GRIDTRANSFER_H_JCMHPZQU + +#include "BlockProcessor_MPI.h" +#include "Cubism/BlockInfo.h" + +#include +#include + +using namespace cubism; + +struct grid_smoother +{ + StencilInfo stencil; + int stencil_start[3]; + int stencil_end[3]; + + grid_smoother() : stencil(-1, -1, -1, 2, 2, 2, false, 1, 0) + { + stencil_start[0] = stencil_start[1] = stencil_start[2] = -1; + stencil_end[0] = stencil_end[1] = stencil_end[2] = 2; + } + + grid_smoother(const grid_smoother &c) + : stencil(-1, -1, -1, 2, 2, 2, false, 1, 0) + { + stencil_start[0] = stencil_start[1] = stencil_start[2] = -1; + stencil_end[0] = stencil_end[1] = stencil_end[2] = 2; + } + + template + inline void operator()(TLab& lab, const BlockInfo& info, TBlock& o) const + { + typedef typename TBlock::ElementType TElement; + for(int iz=0; iz +struct grid_transfer +{ + // define variables + + // number of neighboring (coarse) cells that are + // required to compute the current values at a fine cell + // (corresponds to stencil when discrtizing pdes) + StencilInfo stencil; + int stencil_start[3], stencil_end[3]; + + // point to the fine mesh which should be filled + typedef typename TGrid::BlockType TBlock; + TGrid& fine_grid; + + // number of coarse blocks in each spatial dimension + const int cNX; + const int cNY; + const int cNZ; + + // number of fine blocks in each spatial dimension + const int fNX; + const int fNY; + const int fNZ; + + // contructor + grid_transfer(TGrid &fgrid, + const int cnx, + const int cny, + const int cnz, + const int fnx, + const int fny, + const int fnz) + : stencil(-2, -2, -2, 3, 3, 3, true, 1, 0), fine_grid(fgrid), cNX(cnx), + cNY(cny), cNZ(cnz), fNX(fnx), fNY(fny), fNZ(fnz) + { + stencil_start[0] = stencil_start[1] = stencil_start[2] = -2; + stencil_end[0] = stencil_end[1] = stencil_end[2] = 3; + } + + // copy constructor + grid_transfer(const grid_transfer >) + : stencil(-2, -2, -2, 3, 3, 3, true, 1, 0), fine_grid(gt.fine_grid), + cNX(gt.cNX), cNY(gt.cNY), cNZ(gt.cNZ), fNX(gt.fNX), fNY(gt.fNY), + fNZ(gt.fNZ) + { + stencil_start[0] = stencil_start[1] = stencil_start[2] = -2; + stencil_end[0] = stencil_end[1] = stencil_end[2] = 3; + } + + // function to tranfer x,y,z-index space into block id for fine grid + int calc_ID_fine(const int ix, const int iy, const int iz) + { + return ix + fNX*(iy + fNY*iz); + } + + // function to calculate offset for index space of cell for coarse grid + std::vector calc_idx_offset(const int ix, const int iy, const int iz) + { + std::vector idx_off(3,0); + idx_off[0] = (ix%2) * TBlock::sizeX / 2; + idx_off[1] = (iy%2) * TBlock::sizeY / 2; + idx_off[2] = (iz%2) * TBlock::sizeZ / 2; + return idx_off; + } + + // references: + // - A. Harten, Multiresolution algorithms for the numerical solution of hyperbolic conservation laws, + // Comm. Pur Appl. Math. 48 (1995) 1305-1342. + // - B. L. Bihari & A. Harten, Multiresolution schemes for the numerical solution of of 2-d conservation laws I. + // SIAM J. Sci. Comput. 18 (1997) 315-354. + // - O. Roussel et al., A conservative fully adaptive multiresolution algorithm for parabolic PDEs, + // Journal Comput. Phys. 188 (2003) 493-523. + // Be careful with the signs given in those references!. + template + ReturnType _interpolate_HARTEN(const int f_cell_ix, // first/x index of fine cell + const int f_cell_iy, // second/y index of fine cell + const int f_cell_iz, // third/z index of fine cell + const int c_cell_ix_off, // first/x index offset of coarse cell + const int c_cell_iy_off, // second/y index offset of coarse cell + const int c_cell_iz_off, // third/z index offset of coarse cell + LabType &lab, // the lab corresponding to the coarse cell + const int s, // stencil width related to the order r as r=2s+1 + const Real gamma[2]) // weights + { + // to be filled + ReturnType out; + + // index to distinguish first and second fine cell in coarse cell in one spatial direction + const int n = f_cell_ix%2; + const int p = f_cell_iy%2; + const int q = f_cell_iz%2; + + // get coarse cell indices containing the current fine cell + const int c_cell_ix = c_cell_ix_off + (f_cell_ix - n)/2; + const int c_cell_iy = c_cell_iy_off + (f_cell_iy - p)/2; + const int c_cell_iz = c_cell_iz_off + (f_cell_iz - q)/2; + + // compute all terms + ReturnType Qsx; + Qsx.clear(); + ReturnType Qsy; + Qsy.clear(); + ReturnType Qsz; + Qsz.clear(); + ReturnType Qsxy; + Qsxy.clear(); + ReturnType Qsxz; + Qsxz.clear(); + ReturnType Qsyz; + Qsyz.clear(); + ReturnType Qsxyz; + Qsxyz.clear(); + + for (int rr=1; rr<=s; rr++) + { + Qsx = Qsx + gamma[rr-1] * (lab(c_cell_ix+rr, c_cell_iy, c_cell_iz) - lab(c_cell_ix-rr, c_cell_iy, c_cell_iz)); + Qsy = Qsy + gamma[rr-1] * (lab(c_cell_ix, c_cell_iy+rr, c_cell_iz) - lab(c_cell_ix, c_cell_iy-rr, c_cell_iz)); + Qsz = Qsz + gamma[rr-1] * (lab(c_cell_ix, c_cell_iy, c_cell_iz+rr) - lab(c_cell_ix, c_cell_iy, c_cell_iz-rr)); + + for (int ss=1; ss<=s; ss++) + { + Qsxy = Qsxy + gamma[rr-1] * gamma[ss-1] * (lab(c_cell_ix+rr, c_cell_iy+ss, c_cell_iz) - lab(c_cell_ix+rr, c_cell_iy-ss, c_cell_iz) + - lab(c_cell_ix-rr, c_cell_iy+ss, c_cell_iz) + lab(c_cell_ix-rr, c_cell_iy-ss, c_cell_iz)); + Qsxz = Qsxz + gamma[rr-1] * gamma[ss-1] * (lab(c_cell_ix+rr, c_cell_iy, c_cell_iz+ss) - lab(c_cell_ix+rr, c_cell_iy, c_cell_iz-ss) + - lab(c_cell_ix-rr, c_cell_iy, c_cell_iz+ss) + lab(c_cell_ix-rr, c_cell_iy, c_cell_iz-ss)); + Qsyz = Qsyz + gamma[rr-1] * gamma[ss-1] * (lab(c_cell_ix, c_cell_iy+rr, c_cell_iz+ss) - lab(c_cell_ix, c_cell_iy+rr, c_cell_iz-ss) + - lab(c_cell_ix, c_cell_iy-rr, c_cell_iz+ss) + lab(c_cell_ix, c_cell_iy-rr, c_cell_iz-ss)); + + for (int tt=1; tt<=s; tt++) + Qsxyz = Qsxyz + gamma[rr-1] * gamma[ss-1] * gamma[tt-1] * (lab(c_cell_ix+rr, c_cell_iy+ss, c_cell_iz+tt) - lab(c_cell_ix+rr, c_cell_iy+ss, c_cell_iz-tt) + - lab(c_cell_ix+rr, c_cell_iy-ss, c_cell_iz+tt) - lab(c_cell_ix-rr, c_cell_iy+ss, c_cell_iz+tt) + + lab(c_cell_ix+rr, c_cell_iy-ss, c_cell_iz-tt) + lab(c_cell_ix-rr, c_cell_iy+ss, c_cell_iz-tt) + + lab(c_cell_ix-rr, c_cell_iy-ss, c_cell_iz+tt) - lab(c_cell_ix-rr, c_cell_iy-ss, c_cell_iz-tt)); + } + } + + // sum up + out = lab(c_cell_ix, c_cell_iy, c_cell_iz) + + static_cast(pow(-1.0,(double) n)) * Qsx + static_cast(pow(-1.0,(double) p)) * Qsy + static_cast(pow(-1.0,(double) q)) * Qsz + + static_cast(pow(-1.0,(double) (n+p))) * Qsxy + static_cast(pow(-1.0,(double) (n+q))) * Qsxz + static_cast(pow(-1.0,(double) (p+q))) * Qsyz + + static_cast(pow(-1.0,(double) (n+p+q))) * Qsxyz; + + return out; + } + + + template + void operator()(LabType& lab, const BlockInfo& info, BlockType& o) + { + // get id of current coarse grid block + const long long cID = info.blockID; + + // get coarse block in x,y,z-index space + const int cbiz = (int) floor(((double) cID) / ((double) (cNX * cNY))); + const int tmp = cID % (cNX * cNY); + const int cbiy = (int) floor(((double) tmp) / ((double) (cNX))); + const int cbix = tmp % cNX; + + // get first fine block in x,y,z-index space + const int fbix = 2 * cbix; + const int fbiy = 2 * cbiy; + const int fbiz = 2 * cbiz; + // get vector of all fine blocks in current coarse block + std::vector< std::vector > fidx(8); + std::vector vtmp(3); + vtmp[0] = fbix; vtmp[1] = fbiy; vtmp[2] = fbiz; + fidx[0] = vtmp; + vtmp[0] = fbix+1; vtmp[1] = fbiy; vtmp[2] = fbiz; + fidx[1] = vtmp; + vtmp[0] = fbix; vtmp[1] = fbiy+1; vtmp[2] = fbiz; + fidx[2] = vtmp; + vtmp[0] = fbix; vtmp[1] = fbiy; vtmp[2] = fbiz+1; + fidx[3] = vtmp; + vtmp[0] = fbix+1; vtmp[1] = fbiy+1; vtmp[2] = fbiz; + fidx[4] = vtmp; + vtmp[0] = fbix+1; vtmp[1] = fbiy; vtmp[2] = fbiz+1; + fidx[5] = vtmp; + vtmp[0] = fbix; vtmp[1] = fbiy+1; vtmp[2] = fbiz+1; + fidx[6] = vtmp; + vtmp[0] = fbix+1; vtmp[1] = fbiy+1; vtmp[2] = fbiz+1; + fidx[7] = vtmp; + + // get all blocks of fine grid + std::vector f_info = fine_grid.getBlocksInfo(); + + typedef typename BlockType::ElementType TElement; + + // loop all blocks of fine grid within the current block of the coarse grid + for(int i=0; i<8; i++) + { + // get fine-grid block via its id calculated from indices + const int fID = calc_ID_fine(fidx[i][0], fidx[i][1], fidx[i][2]); + BlockInfo cf_info = f_info[fID]; + BlockType& b = *(BlockType*)cf_info.ptrBlock; + + // calculate starting index for coarse cells comprising current fine grid block + std::vector c_idx_off = calc_idx_offset(fidx[i][0], fidx[i][1], fidx[i][2]); + + // loop all cells of this fine grid block + for (int iz=0; iz(ix, iy, iz, c_idx_off[0], c_idx_off[1], c_idx_off[2], lab, 1, gamma); +#else + // 5th order accurate interpolation + const Real gamma1 = -22.0/128.0; + const Real gamma2 = 3.0/128.0; + const Real gamma[2] = {gamma1, gamma2}; + const TElement actcell = _interpolate_HARTEN(ix, iy, iz, c_idx_off[0], c_idx_off[1], c_idx_off[2], lab, 2, gamma); +#endif + b(ix,iy,iz) = actcell; + } + } + + return; + } + +}; + +template +void interpolate_from_coarse(const TGridIn &coarse_grid, + TGridOut &fine_grid, + const int in_bpdx, + const int in_bpdy, + const int in_bpdz, + const int out_bpdx, + const int out_bpdy, + const int out_bpdz, + const size_t smooth_iter = 0, + const bool verbose = true) +{ + // construct (previous) coarser grid + // we assume that the resolution increases by a factor of 2 + assert(out_bpdx%2 == 0 && out_bpdz%2 == 0 && out_bpdz%2 == 0); + assert(2*in_bpdx == out_bpdx); + assert(2*in_bpdy == out_bpdy); + assert(2*in_bpdz == out_bpdz); + + // do the interpolation from the coarse to the fine mesh + // we fill the fine grid on each rank + // therefore, we have to set up a halo/ghost cell layer for the + // coarse grid (hence, the coarse grid goes into the lab) + grid_transfer coarse_to_fine(fine_grid, + in_bpdx, in_bpdy, in_bpdz, + out_bpdx, out_bpdy, out_bpdz); + + // This requires the that coarse block on THIS process contains the fine + // blocks + // + // o------x------o + // | | | + // | | | + // x------x------x + // | | | + // | | | + // o------x------o + // Coarse block corners "o", fine blocks corners "x", both must be in the + // memory space of this process + process(coarse_to_fine, const_cast(coarse_grid), 0, 0); + + // smooth grid + for (size_t i = 0; i < smooth_iter; ++i) { + if (verbose) + std::cout << "smoothing grid: iteration " << i+1 << std::endl; + grid_smoother smoother; + process(smoother, fine_grid, 0, 0); + } +} + +#endif /* MPI_GRIDTRANSFER_H_JCMHPZQU */ diff --git a/tools/gridmanip/Prolongation/ProlongHarten.h b/tools/gridmanip/Prolongation/ProlongHarten.h new file mode 100644 index 0000000000000000000000000000000000000000..996632b358827d5abb1110715c857d8b860c635e --- /dev/null +++ b/tools/gridmanip/Prolongation/ProlongHarten.h @@ -0,0 +1,61 @@ +// File : ProlongHarten.h +// Created : Tue Jul 11 2017 03:01:34 PM (+0200) +// Author : Fabian Wermelinger +// Description: Prolongation based on Harten (MPI_GridTransfer.h) +// Copyright 2017 ETH Zurich. All Rights Reserved. +#ifndef PROLONGHARTEN_H_AFYYD4ZZ +#define PROLONGHARTEN_H_AFYYD4ZZ + +#include "GridOperator.h" +#include "Prolongation/MPI_GridTransfer.h" + +#include +#include + +using namespace cubism; + +template +class ProlongHarten : public GridOperator +{ +public: + ProlongHarten(ArgumentParser &p) + : GridOperator(p) + { + } + + ~ProlongHarten() = default; + + void operator()(const TGridIn &grid_in, + TGridOut &grid_out, + const bool verbose) override + { + // 0.) checks + assert(TGridIn::BlockType::sizeX == TGridOut::BlockType::sizeX); + assert(TGridIn::BlockType::sizeY == TGridOut::BlockType::sizeY); + assert(TGridIn::BlockType::sizeZ == TGridOut::BlockType::sizeZ); + const int NX_in = grid_in.getResidentBlocksPerDimension(0); + const int NY_in = grid_in.getResidentBlocksPerDimension(1); + const int NZ_in = grid_in.getResidentBlocksPerDimension(2); + const int NX_out = grid_out.getResidentBlocksPerDimension(0); + const int NY_out = grid_out.getResidentBlocksPerDimension(1); + const int NZ_out = grid_out.getResidentBlocksPerDimension(2); + assert(NX_in <= NX_out); + assert(NY_in <= NY_out); + assert(NZ_in <= NZ_out); + + const size_t smooth_iter = this->m_parser("smooth_iter").asInt(0); + + interpolate_from_coarse(grid_in, + grid_out, + NX_in, + NY_in, + NZ_in, + NX_out, + NY_out, + NZ_out, + smooth_iter, + verbose); + } +}; + +#endif /* PROLONGHARTEN_H_AFYYD4ZZ */ diff --git a/tools/gridmanip/README.md b/tools/gridmanip/README.md new file mode 100644 index 0000000000000000000000000000000000000000..058be0b7ba96fea6c6e54cfad49cb907b6692631 --- /dev/null +++ b/tools/gridmanip/README.md @@ -0,0 +1,22 @@ +# Grid manipulation tool + +The tool uses the I/O routines provided by the Cubism library to perform +various grid operations. If necessary, additional I/O routines can be added +to support non-native data formats. + +## Prolongation + +Computes a refined grid using the interpolation scheme suggested by Harten et +al. (1995 and 1997), see also the file `Prolongation/MPI_GridTransfer.h`. +Refinement steps that are a power of two of the input grid are supported at the +moment. + +## Restriction + +Computes a coarsening of the provided input grid by averaging adjacent cells +into a coarser cell. Coarsening by a power of two is supported at the moment. + +## Smoothing + +This operation performs a number of smoothing steps on based on the data of the +input grid. A 7-point or 27-point stencil is supported at the moment. diff --git a/tools/gridmanip/Restriction/RestrictBlockAverage.h b/tools/gridmanip/Restriction/RestrictBlockAverage.h new file mode 100644 index 0000000000000000000000000000000000000000..4aabc3bf9254995e58b60e7753f3cc7e5b7788c0 --- /dev/null +++ b/tools/gridmanip/Restriction/RestrictBlockAverage.h @@ -0,0 +1,115 @@ +// File : RestrictBlockAverage.h +// Created : Mon Jul 10 2017 12:41:30 PM (+0200) +// Author : Fabian Wermelinger +// Description: Restriction operator by simple averaging +// Copyright 2017 ETH Zurich. All Rights Reserved. +#ifndef RESTRICTBLOCKAVERAGE_H_EJPTMN1I +#define RESTRICTBLOCKAVERAGE_H_EJPTMN1I + +#include "Cubism/BlockInfo.h" +#include "GridOperator.h" + +#include +#include + +using namespace cubism; + +template +class RestrictBlockAverage : public GridOperator +{ +public: + RestrictBlockAverage(ArgumentParser &p) + : GridOperator(p) + { + } + ~RestrictBlockAverage() = default; + + void operator()(const TGridIn &grid_in, + TGridOut &grid_out, + const bool verbose) override + { + // 0.) checks + typedef typename TGridIn::BlockType TBlockIn; + typedef typename TGridOut::BlockType TBlockOut; + assert(TBlockIn::sizeX == TBlockOut::sizeX); + assert(TBlockIn::sizeY == TBlockOut::sizeY); + assert(TBlockIn::sizeZ == TBlockOut::sizeZ); + const int NX_in = grid_in.getResidentBlocksPerDimension(0); + const int NY_in = grid_in.getResidentBlocksPerDimension(1); + const int NZ_in = grid_in.getResidentBlocksPerDimension(2); + const int NX_out = grid_out.getResidentBlocksPerDimension(0); + const int NY_out = grid_out.getResidentBlocksPerDimension(1); + const int NZ_out = grid_out.getResidentBlocksPerDimension(2); + assert(NX_in >= NX_out); // our goal is to coarsen + assert(NY_in >= NY_out); // our goal is to coarsen + assert(NZ_in >= NZ_out); // our goal is to coarsen + + // 1.) assign fine to coarse map + assert(NX_in % NX_out == 0); + assert(NY_in % NY_out == 0); + assert(NZ_in % NZ_out == 0); + const int cpdx = NX_in / NX_out; + const int cpdy = NY_in / NY_out; + const int cpdz = NZ_in / NZ_out; + assert(cpdx > 0 && cpdx <= TBlockOut::sizeX); + assert(cpdy > 0 && cpdy <= TBlockOut::sizeY); + assert(cpdz > 0 && cpdz <= TBlockOut::sizeZ); + + assert(TBlockOut::sizeX % cpdx == 0); + assert(TBlockOut::sizeY % cpdy == 0); + assert(TBlockOut::sizeZ % cpdz == 0); + const int bstrideX = TBlockOut::sizeX / cpdx; + const int bstrideY = TBlockOut::sizeY / cpdy; + const int bstrideZ = TBlockOut::sizeZ / cpdz; + + // 2.) average + std::vector info_in = grid_in.getResidentBlocksInfo(); + std::vector info_out = grid_out.getResidentBlocksInfo(); + +#pragma omp parallel for + for (size_t i = 0; i < info_out.size(); i++) { + BlockInfo infoout = info_out[i]; + TBlockOut &bout = *(TBlockOut *)infoout.ptrBlock; + bout.clear(); // zero data + } + + for (size_t i = 0; i < info_in.size(); i++) { + // src + BlockInfo infoin = info_in[i]; + TBlockIn &bin = *(TBlockIn *)infoin.ptrBlock; + + // dst + const int bix = infoin.index[0] / cpdx; + const int biy = infoin.index[1] / cpdy; + const int biz = infoin.index[2] / cpdz; + const int cox = bstrideX * (infoin.index[0] % cpdx); + const int coy = bstrideY * (infoin.index[1] % cpdy); + const int coz = bstrideZ * (infoin.index[2] % cpdz); + BlockInfo infoout = info_out[bix + NX_out * (biy + NY_out * biz)]; + TBlockOut &bout = *(TBlockOut *)infoout.ptrBlock; + + for (int iz = 0; iz < TBlockIn::sizeZ; iz++) + for (int iy = 0; iy < TBlockIn::sizeY; iy++) + for (int ix = 0; ix < TBlockIn::sizeX; ix++) { + const int cix = ix / cpdx + cox; + const int ciy = iy / cpdy + coy; + const int ciz = iz / cpdz + coz; + bout(cix, ciy, ciz) = + bout(cix, ciy, ciz) + bin(ix, iy, iz); + } + } + + const Real factor = 1.0 / (cpdx * cpdy * cpdz); +#pragma omp parallel for + for (size_t i = 0; i < info_out.size(); i++) { + BlockInfo infoout = info_out[i]; + TBlockOut &bout = *(TBlockOut *)infoout.ptrBlock; + for (int iz = 0; iz < TBlockOut::sizeZ; iz++) + for (int iy = 0; iy < TBlockOut::sizeY; iy++) + for (int ix = 0; ix < TBlockOut::sizeX; ix++) + bout(ix, iy, iz) = factor * bout(ix, iy, iz); + } + } +}; + +#endif /* RESTRICTBLOCKAVERAGE_H_EJPTMN1I */ diff --git a/tools/gridmanip/Smoother/Smoother.h b/tools/gridmanip/Smoother/Smoother.h new file mode 100644 index 0000000000000000000000000000000000000000..486b5f8176bec0fe7d6cdea928ff4e08f97ddc24 --- /dev/null +++ b/tools/gridmanip/Smoother/Smoother.h @@ -0,0 +1,85 @@ +// File : Smoother.h +// Created : Sun Oct 29 2017 12:36:21 PM (+0100) +// Author : Fabian Wermelinger +// Description: Smooth data +// Copyright 2017 ETH Zurich. All Rights Reserved. +#ifndef SMOOTHER_H_2PBGUG6D +#define SMOOTHER_H_2PBGUG6D + +#include "Cubism/BlockInfo.h" +#include "GridOperator.h" +#include "Prolongation/MPI_GridTransfer.h" + +#include +#include + +using namespace cubism; + +template +class Smoother : public GridOperator +{ +public: + Smoother(ArgumentParser &p) : GridOperator(p) + { + } + + ~Smoother() = default; + + void operator()(const TGridIn &grid_in, + TGridOut &grid_out, + const bool verbose) override + { + // 0.) checks + typedef typename TGridIn::BlockType TBlockIn; + typedef typename TGridOut::BlockType TBlockOut; + assert(TBlockIn::sizeX == TBlockOut::sizeX); + assert(TBlockIn::sizeY == TBlockOut::sizeY); + assert(TBlockIn::sizeZ == TBlockOut::sizeZ); + assert(grid_in.getResidentBlocksPerDimension(0) == + grid_out.getResidentBlocksPerDimension(0)); + assert(grid_in.getResidentBlocksPerDimension(1) == + grid_out.getResidentBlocksPerDimension(1)); + assert(grid_in.getResidentBlocksPerDimension(2) == + grid_out.getResidentBlocksPerDimension(2)); + + const size_t smooth_iter = this->m_parser("smooth_iter").asInt(0); + + // copy over + std::vector info_in = grid_in.getResidentBlocksInfo(); + std::vector info_out = grid_out.getResidentBlocksInfo(); + assert(info_in.size() == info_out.size()); + +#pragma omp parallel for + for (size_t i = 0; i < info_out.size(); i++) { + BlockInfo infoout = info_out[i]; + TBlockOut &bout = *(TBlockOut *)infoout.ptrBlock; + bout.clear(); // zero data + } + +#pragma omp parallel for + for (size_t i = 0; i < info_in.size(); i++) { + // src + BlockInfo infoin = info_in[i]; + TBlockIn &bin = *(TBlockIn *)infoin.ptrBlock; + + // dst + BlockInfo infoout = info_out[i]; + TBlockOut &bout = *(TBlockOut *)infoout.ptrBlock; + + for (int iz = 0; iz < TBlockIn::sizeZ; iz++) + for (int iy = 0; iy < TBlockIn::sizeY; iy++) + for (int ix = 0; ix < TBlockIn::sizeX; ix++) + bout(ix, iy, iz) = bin(ix, iy, iz); + } + + // smooth out grid + for (size_t i = 0; i < smooth_iter; ++i) { + if (verbose) + std::cout << "smoothing grid: iteration " << i + 1 << std::endl; + grid_smoother smoother; + process(smoother, grid_out, 0, 0); + } + } +}; + +#endif /* SMOOTHER_H_2PBGUG6D */ diff --git a/tools/gridmanip/Types.h b/tools/gridmanip/Types.h new file mode 100644 index 0000000000000000000000000000000000000000..54538607c5c3a04c1239743ed480a220168df604 --- /dev/null +++ b/tools/gridmanip/Types.h @@ -0,0 +1,389 @@ +// File : Types.h +// Created : Sat May 04 2019 09:41:16 PM (+0200) +// Author : Fabian Wermelinger +// Description: Type declarations for grid and boundary processing +// Copyright 2019 ETH Zurich. All Rights Reserved. +#ifndef TYPES_H_AESHUI1N +#define TYPES_H_AESHUI1N + +#ifdef _FLOAT_PRECISION_ +typedef float Real; +typedef float DumpReal; +#else +typedef double Real; +typedef double DumpReal; +#endif + +// scalar field pass-through streamer. assumes that the element type TElement +// is a struct with field member .scalar. +struct StreamerScalar { + static const std::string NAME; + static const std::string EXT; + static const int NCHANNELS = 1; + static const int CLASS = 0; + + template + static inline void operate(const TBlock &b, + const int ix, + const int iy, + const int iz, + TReal output[NCHANNELS]) + { + typedef typename TBlock::ElementType TElement; + const TElement &el = b(ix, iy, iz); + output[0] = el.scalar; + } + + template + static inline void operate(TBlock &b, + const TReal input[NCHANNELS], + const int ix, + const int iy, + const int iz) + { + typedef typename TBlock::ElementType TElement; + TElement &el = b(ix, iy, iz); + el.scalar = input[0]; + } + + static const char *getAttributeName() { return "Scalar"; } +}; + +const std::string StreamerScalar::NAME = "Scalar"; +const std::string StreamerScalar::EXT = "-scalar"; + +struct GridElement { + using RealType = Real; + + RealType scalar; + + void clear() { scalar = 0.0; } + + GridElement &operator=(const GridElement &gp) + { + this->scalar = gp.scalar; + return *this; + } + + GridElement &operator+=(const GridElement &rhs) + { + scalar += rhs.scalar; + return *this; + } + + GridElement &operator-=(const GridElement &rhs) + { + scalar -= rhs.scalar; + return *this; + } + + GridElement &operator*=(const RealType c) + { + scalar *= c; + return *this; + } + + friend GridElement operator+(GridElement lhs, const GridElement &rhs) + { + return (lhs += rhs); + } + + friend GridElement operator-(GridElement lhs, const GridElement &rhs) + { + return (lhs -= rhs); + } + + friend GridElement operator*(GridElement lhs, const RealType c) + { + return (lhs *= c); + } + + friend GridElement operator*(const RealType c, GridElement rhs) + { + return (rhs *= c); + } +}; + +struct GridBlock { + static const int sizeX = _BLOCKSIZE_; + static const int sizeY = _BLOCKSIZE_; + static const int sizeZ = _BLOCKSIZE_; + + static const int gptfloats = + sizeof(GridElement) / sizeof(GridElement::RealType); + + using RealType = typename GridElement::RealType; + using ElementType = GridElement; + using element_type = GridElement; + + GridElement __attribute__((__aligned__(CUBISM_ALIGNMENT))) + data[_BLOCKSIZE_][_BLOCKSIZE_][_BLOCKSIZE_]; + + RealType __attribute__((__aligned__(CUBISM_ALIGNMENT))) + tmp[_BLOCKSIZE_][_BLOCKSIZE_][_BLOCKSIZE_][gptfloats]; + + void clear_data() + { + const int N = sizeX * sizeY * sizeZ; + GridElement *const e = &data[0][0][0]; + for (int i = 0; i < N; ++i) { + e[i].clear(); + } + } + + void clear_tmp() + { + const int N = sizeX * sizeY * sizeZ * gptfloats; + + RealType *const e = &tmp[0][0][0][0]; + for (int i = 0; i < N; ++i) { + e[i] = 0; + } + } + + void clear() + { + clear_data(); + clear_tmp(); + } + + inline GridElement &operator()(int ix, int iy = 0, int iz = 0) + { + assert(ix >= 0 && ix < sizeX); + assert(iy >= 0 && iy < sizeY); + assert(iz >= 0 && iz < sizeZ); + return data[iz][iy][ix]; + } + + inline const GridElement &operator()(int ix, int iy = 0, int iz = 0) const + { + assert(ix >= 0 && ix < sizeX); + assert(iy >= 0 && iy < sizeY); + assert(iz >= 0 && iz < sizeZ); + return data[iz][iy][ix]; + } +}; + +template class Alloc = std::allocator> +class ExtrapolatingBoundaryTensorial : public cubism::BlockLab +{ + using ElementTypeBlock = typename TBlock::ElementType; + +public: + std::string name() const override + { + return "ExtrapolatingBoundaryTensorial"; + } + bool is_xperiodic() override { return false; } + bool is_yperiodic() override { return false; } + bool is_zperiodic() override { return false; } + + ExtrapolatingBoundaryTensorial() = default; + + void _apply_bc(const cubism::BlockInfo &info, const Real t = 0) override + { + if (info.index[0] == 0) + this->template extrapolate_<0, 0>(); + if (info.index[0] == this->NX - 1) + this->template extrapolate_<0, 1>(); + if (info.index[1] == 0) + this->template extrapolate_<1, 0>(); + if (info.index[1] == this->NY - 1) + this->template extrapolate_<1, 1>(); + if (info.index[2] == 0) + this->template extrapolate_<2, 0>(); + if (info.index[2] == this->NZ - 1) + this->template extrapolate_<2, 1>(); + } + +private: + int start_[3], end_[3]; + + template + void setupStartEnd_() + { + start_[0] = dir == 0 + ? (side == 0 ? this->m_stencilStart[0] : TBlock::sizeX) + : 0; + start_[1] = dir == 1 + ? (side == 0 ? this->m_stencilStart[1] : TBlock::sizeY) + : 0; + start_[2] = dir == 2 + ? (side == 0 ? this->m_stencilStart[2] : TBlock::sizeZ) + : 0; + + end_[0] = + dir == 0 + ? (side == 0 ? 0 : TBlock::sizeX + this->m_stencilEnd[0] - 1) + : TBlock::sizeX; + end_[1] = + dir == 1 + ? (side == 0 ? 0 : TBlock::sizeY + this->m_stencilEnd[1] - 1) + : TBlock::sizeY; + end_[2] = + dir == 2 + ? (side == 0 ? 0 : TBlock::sizeZ + this->m_stencilEnd[2] - 1) + : TBlock::sizeZ; + } + + ElementTypeBlock &accessCacheBlock_(int ix, int iy, int iz) + { + return this->m_cacheBlock->Access(ix - this->m_stencilStart[0], + iy - this->m_stencilStart[1], + iz - this->m_stencilStart[2]); + } + + template + void extrapolate_() + { + this->template setupStartEnd_(); + + // faces + for (int iz = start_[2]; iz < end_[2]; iz++) + for (int iy = start_[1]; iy < end_[1]; iy++) + for (int ix = start_[0]; ix < end_[0]; ix++) { + accessCacheBlock_(ix, iy, iz) = accessCacheBlock_( + dir == 0 ? (side == 0 ? 0 : TBlock::sizeX - 1) : ix, + dir == 1 ? (side == 0 ? 0 : TBlock::sizeY - 1) : iy, + dir == 2 ? (side == 0 ? 0 : TBlock::sizeZ - 1) : iz); + } + + // edges and corners + { + int s_[3], e_[3]; + const int bsize[3] = {TBlock::sizeX, TBlock::sizeY, TBlock::sizeZ}; + + s_[dir] = + this->m_stencilStart[dir] * (1 - side) + bsize[dir] * side; + e_[dir] = (bsize[dir] - 1 + this->m_stencilEnd[dir]) * side; + + const int d1 = (dir + 1) % 3; + const int d2 = (dir + 2) % 3; + + for (int b = 0; b < 2; ++b) + for (int a = 0; a < 2; ++a) { + s_[d1] = this->m_stencilStart[d1] + + a * b * (bsize[d1] - this->m_stencilStart[d1]); + s_[d2] = + this->m_stencilStart[d2] + + (a - a * b) * (bsize[d2] - this->m_stencilStart[d2]); + + e_[d1] = (1 - b + a * b) * + (bsize[d1] - 1 + this->m_stencilEnd[d1]); + e_[d2] = (a + b - a * b) * + (bsize[d2] - 1 + this->m_stencilEnd[d2]); + + for (int iz = s_[2]; iz < e_[2]; iz++) + for (int iy = s_[1]; iy < e_[1]; iy++) + for (int ix = s_[0]; ix < e_[0]; ix++) { + accessCacheBlock_(ix, iy, iz) = + dir == 0 + ? accessCacheBlock_( + side * (bsize[0] - 1), iy, iz) + : (dir == 1 + ? accessCacheBlock_( + ix, + side * (bsize[1] - 1), + iz) + : accessCacheBlock_( + ix, + iy, + side * (bsize[2] - 1))); + } + } + } + } +}; + +enum class IOType { In = 0, Out }; + +template +class InputOutputFactory +{ +public: + InputOutputFactory(const bool isroot, cubism::ArgumentParser &p) + : m_isroot(isroot), m_parser(p), finput(nullptr), foutput(nullptr) + { + if (Direction == IOType::In) { + setInputStreamer_(); + } else if (Direction == IOType::Out) { + setOutputStreamer_(); + } else { + if (m_isroot) + std::cerr << "ERROR: Unknown Direction type" << std::endl; + abort(); + } + } + + using InStreamer = void (*)(TGrid &, + const std::string &, + const std::string &); + using OutStreamer = void (*)(const TGrid &, + const int, + const Real, + const std::string &, + const std::string &, + const bool); + + void operator()(TGrid &grid, + const std::string &fname, + const std::string &path = ".") + { + if (Direction == IOType::In) { + finput(grid, fname, path); + } else { + foutput(grid, 0, 0.0, fname, path, true); + } + } + +private: + const bool m_isroot; + cubism::ArgumentParser &m_parser; + + InStreamer finput; + OutStreamer foutput; + + void setOutputStreamer_() + { + const std::string save_format = + m_parser("save_format").asString("cubism_h5"); + + if (save_format == "cubism_h5") { + foutput = &cubism::DumpHDF5_MPI; + } else { + if (m_isroot) + std::cerr << "ERROR: No suitable save format chosen" + << std::endl; + abort(); + } + + if (foutput == nullptr) { + if (m_isroot) + std::cerr << "ERROR: NULL kernel" << std::endl; + abort(); + } + } + + void setInputStreamer_() + { + const std::string read_format = + m_parser("read_format").asString("cubism_h5"); + + if (read_format == "cubism_h5") { + finput = &cubism::ReadHDF5_MPI; + } else { + if (m_isroot) + std::cerr << "ERROR: No suitable input format chosen" + << std::endl; + abort(); + } + + if (finput == nullptr) { + if (m_isroot) + std::cerr << "ERROR: NULL kernel" << std::endl; + abort(); + } + } +}; + +#endif /* TYPES_H_AESHUI1N */ diff --git a/tools/gridmanip/main.cpp b/tools/gridmanip/main.cpp new file mode 100644 index 0000000000000000000000000000000000000000..6e3fc223735cb2cf0dcec083a95fbb993337a4db --- /dev/null +++ b/tools/gridmanip/main.cpp @@ -0,0 +1,136 @@ +// File : main.cpp +// Created : Mon Jul 10 2017 12:18:47 PM (+0200) +// Author : Fabian Wermelinger +// Description: Tool for restriction / prolongation of grid +// Copyright 2017 ETH Zurich. All Rights Reserved. +#ifndef CUBISM_USE_HDF +#pragma error "HDF5 library is required for build" +#endif /* CUBISM_USE_HDF */ + +#include "Cubism/ArgumentParser.h" +#include "Cubism/BlockInfo.h" +#include "Cubism/BlockLab.h" +#include "Cubism/BlockLabMPI.h" +#include "Cubism/Grid.h" +#include "Cubism/GridMPI.h" +#include "Cubism/HDF5Dumper_MPI.h" + +// grid operators +#include "GridOperator.h" +#include "Prolongation/ProlongHarten.h" +#include "Restriction/RestrictBlockAverage.h" +#include "Smoother/Smoother.h" + +// types required to process grids +#include "Types.h" + +#include +#include +#include +#include +#include + +using namespace cubism; + +int main(int argc, char *argv[]) +{ + // Boundary conditions for prolongation (requires tensorial stencil) and + // smoother (conventional 7-point stencil): + // + // periodic (naturally tensorial) + // using LabMPI = BlockLabMPI>; + // + // zeroth-order extrapolating boundary (copy first interior cell into halos) + using LabMPI = BlockLabMPI>; + + int provided; + MPI_Init_thread(&argc, (char ***)&argv, MPI_THREAD_MULTIPLE, &provided); + + int rank; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + const bool isroot = (0 == rank) ? true : false; + + ArgumentParser parser(argc, argv); + if (parser.exist("conf")) + parser.readFile(parser("conf").asString()); + + // setup grids + const Real extent = parser("extent").asDouble(1.0); + const int xpesize = parser("xpesize").asInt(1); + const int ypesize = parser("ypesize").asInt(1); + const int zpesize = parser("zpesize").asInt(1); + + // fine blocks + const int in_bpdx = parser("in_bpdx").asInt(1); + const int in_bpdy = parser("in_bpdy").asInt(1); + const int in_bpdz = parser("in_bpdz").asInt(1); + + // coarse blocks + const int out_bpdx = parser("out_bpdx").asInt(1); + const int out_bpdy = parser("out_bpdy").asInt(1); + const int out_bpdz = parser("out_bpdz").asInt(1); + + using TGridIn = GridMPI>; + using TGridOut = GridMPI>; + + TGridIn *const grid_in = new TGridIn(xpesize, + ypesize, + zpesize, + in_bpdx, + in_bpdy, + in_bpdz, + extent, + MPI_COMM_WORLD); + TGridOut *const grid_out = new TGridOut(xpesize, + ypesize, + zpesize, + out_bpdx, + out_bpdy, + out_bpdz, + extent, + MPI_COMM_WORLD); + + // setup I/O + InputOutputFactory infactory(isroot, parser); + InputOutputFactory outfactory(isroot, parser); + + // get data for input grid + parser.set_strict_mode(); + const std::string infile = parser("in_file").asString(); + parser.unset_strict_mode(); + infactory(*grid_in, infile); + + // perform grid manipulation + GridOperator *gridmanip; + const std::string operator_name = + parser("operator").asString("RestrictBlockAverage"); + if (operator_name == "RestrictBlockAverage") + gridmanip = new RestrictBlockAverage(parser); + else if (operator_name == "ProlongHarten") + gridmanip = new ProlongHarten(parser); + else if (operator_name == "Smoother") + gridmanip = new Smoother(parser); + else { + if (isroot) + std::cerr << "ERROR: Undefined operator '" << operator_name << "'" + << std::endl; + abort(); + } + if (isroot) + parser.print_args(); + (*gridmanip)(*grid_in, *grid_out, isroot); + + // save manipulations + const std::string outfile = parser("out_file").asString("scalar_out"); + outfactory(*grid_out, outfile); + + MPI_Barrier(MPI_COMM_WORLD); + + // clean up + delete grid_in; + delete grid_out; + delete gridmanip; + + MPI_Finalize(); + return 0; +}