Skip to content

GitLab

C
Cubism
Commit ad2e5ac9 authored by fabianw's avatar fabianw
Browse files
Options

Merge branch 'gridmanip'

parents dab6442d 7b9326f6
Show whitespace changes
Inline Side-by-side
Showing with 1482 additions and 0 deletions
tools/gridmanip/.gitignore 0 → 100644
View file @ ad2e5ac9
gridmanip
*.o
*.a
*~
*.swp
*.bak
*.h5
*.xmf
tools/gridmanip/BlockProcessor_MPI.h 0 → 100644
View file @ ad2e5ac9
// 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 <mpi.h>
#include <vector>
using namespace cubism;
template <typename TLab, typename TKernel, typename TGrid>
inline void
process(TKernel rhs, TGrid &grid, const Real t = 0.0, const bool record = false)
{
TKernel myrhs = rhs;
SynchronizerMPI<Real> &Synch = grid.sync(myrhs);
std::vector<BlockInfo> 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 */
tools/gridmanip/Example/Makefile 0 → 100644
View file @ ad2e5ac9
# 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
tools/gridmanip/Example/data/nonuniform/data_000000-5.xmf 0 → 100644
View file @ ad2e5ac9
<?xml version="1.0" ?>
<!DOCTYPE Xdmf SYSTEM "Xdmf.dtd" []>
<Xdmf Version="2.0">
<Domain>
<Grid GridType="Uniform">
<Time Value="0.000000e+00"/>
<Topology TopologyType="3DRectMesh" Dimensions="65 65 65"/>
<Geometry GeometryType="VxVyVz">
<DataItem Name="mesh_vx" Dimensions="65" NumberType="Float" Precision="8" Format="HDF">
data_000000-5.h5:/vx
</DataItem>
<DataItem Name="mesh_vy" Dimensions="65" NumberType="Float" Precision="8" Format="HDF">
data_000000-5.h5:/vy
</DataItem>
<DataItem Name="mesh_vz" Dimensions="65" NumberType="Float" Precision="8" Format="HDF">
data_000000-5.h5:/vz
</DataItem>
</Geometry>
<Attribute Name="data" AttributeType="Scalar" Center="Cell">
<DataItem Dimensions="64 64 64 1" NumberType="Float" Precision="4" Format="HDF">
data_000000-5.h5:/data
</DataItem>
</Attribute>
</Grid>
</Domain>
</Xdmf>
tools/gridmanip/Example/data/uniform/data_000000-5.xmf 0 → 100644
View file @ ad2e5ac9
<?xml version="1.0" ?>
<!DOCTYPE Xdmf SYSTEM "Xdmf.dtd" []>
<Xdmf Version="2.0">
<Domain>
<Grid GridType="Uniform">
<Time Value="0.000000e+00"/>
<Topology TopologyType="3DRectMesh" Dimensions="65 65 65"/>
<Geometry GeometryType="VxVyVz">
<DataItem Name="mesh_vx" Dimensions="65" NumberType="Float" Precision="8" Format="HDF">
data_000000-5.h5:/vx
</DataItem>
<DataItem Name="mesh_vy" Dimensions="65" NumberType="Float" Precision="8" Format="HDF">
data_000000-5.h5:/vy
</DataItem>
<DataItem Name="mesh_vz" Dimensions="65" NumberType="Float" Precision="8" Format="HDF">
data_000000-5.h5:/vz
</DataItem>
</Geometry>
<Attribute Name="data" AttributeType="Scalar" Center="Cell">
<DataItem Dimensions="64 64 64 1" NumberType="Float" Precision="4" Format="HDF">
data_000000-5.h5:/data
</DataItem>
</Attribute>
</Grid>
</Domain>
</Xdmf>
tools/gridmanip/Example/gridmanip 0 → 120000
View file @ ad2e5ac9
../gridmanip
\ No newline at end of file
tools/gridmanip/Example/prolong.conf 0 → 100644
View file @ ad2e5ac9
# 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
tools/gridmanip/Example/restrict.conf 0 → 100644
View file @ ad2e5ac9
# 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
tools/gridmanip/Example/smoother.conf 0 → 100644
View file @ ad2e5ac9
# 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
tools/gridmanip/GridOperator.h 0 → 100644
View file @ ad2e5ac9
// 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 <typename TGridIn, typename TGridOut, typename TBlockLab>
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 */
tools/gridmanip/Makefile 0 → 100644
View file @ ad2e5ac9
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 *~
tools/gridmanip/Prolongation/MPI_GridTransfer.h 0 → 100644
View file @ ad2e5ac9
/* 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 <iostream>
#include <vector>
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<typename TLab, typename TBlock>
inline void operator()(TLab& lab, const BlockInfo& info, TBlock& o) const
{
typedef typename TBlock::ElementType TElement;
for(int iz=0; iz<TBlock::sizeZ; iz++)
for(int iy=0; iy<TBlock::sizeY; iy++)
for(int ix=0; ix<TBlock::sizeX; ix++)
{
// myself
TElement sum = lab(ix,iy,iz);
// neighbor faces
sum = sum + 0.5*(
lab(ix-1,iy,iz) + lab(ix+1,iy,iz) +
lab(ix,iy-1,iz) + lab(ix,iy+1,iz) +
lab(ix,iy,iz-1) + lab(ix,iy,iz+1)
);
// // tensorial
// // neighbor edges
// sum = sum + 0.25*(
// lab(ix,iy-1,iz-1) + lab(ix,iy-1,iz+1) + lab(ix-1,iy-1,iz) + lab(ix+1,iy-1,iz) +
// lab(ix,iy+1,iz-1) + lab(ix,iy+1,iz+1) + lab(ix-1,iy+1,iz) + lab(ix+1,iy+1,iz) +
// lab(ix-1,iy,iz-1) + lab(ix-1,iy,iz+1) + lab(ix+1,iy,iz-1) + lab(ix+1,iy,iz+1)
// );
// // neighbor corners
// sum = sum + 0.125*(
// lab(ix-1,iy-1,iz-1) + lab(ix-1,iy-1,iz+1) + lab(ix-1,iy+1,iz-1) + lab(ix-1,iy+1,iz+1) +
// lab(ix+1,iy-1,iz-1) + lab(ix+1,iy-1,iz+1) + lab(ix+1,iy+1,iz-1) + lab(ix+1,iy+1,iz+1)
// );
// o(ix,iy,iz) = 0.125*sum;
o(ix,iy,iz) = 0.25*sum;
}
}
};
template <typename TGrid>
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 &gt)
: 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<int> calc_idx_offset(const int ix, const int iy, const int iz)
{
std::vector<int> 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<typename LabType, typename ReturnType>
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<Real>(pow(-1.0,(double) n)) * Qsx + static_cast<Real>(pow(-1.0,(double) p)) * Qsy + static_cast<Real>(pow(-1.0,(double) q)) * Qsz
+ static_cast<Real>(pow(-1.0,(double) (n+p))) * Qsxy + static_cast<Real>(pow(-1.0,(double) (n+q))) * Qsxz + static_cast<Real>(pow(-1.0,(double) (p+q))) * Qsyz
+ static_cast<Real>(pow(-1.0,(double) (n+p+q))) * Qsxyz;
return out;
}
template<typename LabType, typename BlockType>
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<int> > fidx(8);
std::vector<int> 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<BlockInfo> 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;