Code indexing in gitaly is broken and leads to code not being visible to the user. We work on the issue with highest priority.

Skip to content
Snippets Groups Projects
Commit 0972cb5b authored by frey_m's avatar frey_m
Browse files

AMR-Test: Remove unnecessary files 

deleted:    ippl/test/AMR/trilinos/AmrTrilinos.cpp
deleted:    ippl/test/AMR/trilinos/AmrTrilinos.h
deleted:    ippl/test/AMR/trilinos/EllipticDomain.cpp
deleted:    ippl/test/AMR/trilinos/EllipticDomain.h
deleted:    ippl/test/AMR/trilinos/IrregularDomain.h
deleted:    ippl/test/AMR/trilinos/testTrilinosGeometry.cpp
parent 309a1498
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
ippl/test/AMR/trilinos/AmrTrilinos.cpp deleted 100644 → 0
+ 0
− 1209
View file @ 309a1498
This diff is collapsed.
ippl/test/AMR/trilinos/AmrTrilinos.h deleted 100644 → 0
+ 0
− 396
View file @ 309a1498
#ifndef AMR_TRILINOS_H
#define AMR_TRILINOS_H
#include <AMReX.H>
#include <AMReX_MultiFab.H>
#include <AMReX_Array.H>
#include <Epetra_MpiComm.h>
#include <Epetra_Map.h>
#include <Epetra_Vector.h>
#include <Epetra_CrsMatrix.h>
#include <Teuchos_RCP.hpp>
#include <Teuchos_ArrayRCP.hpp>
#include <map>
#include <memory>
#include "EllipticDomain.h"
#include "electrostatic_pic_F.H"
#include <AMReX_BCUtil.H>
#include <AMReX_MacBndry.H>
#include <AMReX_FluxRegister.H>
using namespace amrex;
class AmrTrilinos {
public:
typedef std::unique_ptr<MultiFab> AmrField_t;
typedef Array<AmrField_t> container_t;
public:
AmrTrilinos();
void solveWithGeometry(const container_t& rho,
container_t& phi,
container_t& efield,
const Array<Geometry>& geom,
int lbase,
int lfine);
void solve(const container_t& rho,
container_t& phi,
container_t& efield,
const Array<Geometry>& geom,
int lbase,
int lfine,
double scale);
private:
void solveLevel_m(AmrField_t& phi,
const AmrField_t& rho,
const Geometry& geom, int lev);
// single level solve
void linSysSolve_m(AmrField_t& phi, const Geometry& geom);
// void solveGrid_m(
// create map + right-hand side
void build_m(const AmrField_t& rho, const AmrField_t& phi, const Geometry& geom, int lev);
void fillMatrix_m(const Geometry& geom, const AmrField_t& phi, int lev);
int serialize_m(const IntVect& iv) const;
// IntVect deserialize_m(int idx) const;
bool isInside_m(const IntVect& iv) const;
// bool isBoundary_m(const IntVect& iv) const;
// void findBoundaryBoxes_m(const BoxArray& ba);
void copyBack_m(AmrField_t& phi,
const Teuchos::RCP<Epetra_MultiVector>& sol,
const Geometry& geom);
// void grid_m(AmrField_t& rhs, const container_t& rho,
// const Geometry& geom, int lev);
void interpFromCoarseLevel_m(container_t& phi,
const Array<Geometry>& geom,
int lev);
// efield = -grad phi
void getGradient(AmrField_t& phi,
AmrField_t& efield,
const Geometry& geom, int lev);
void getGradient_1(AmrField_t& phi, AmrField_t& efield,
const Geometry& geom, int lev);
void buildLevelMask_m(const BoxArray& grids,
const DistributionMapping& dmap,
Geometry& geom,
int lev)
{
masks_m[lev].reset(new FabArray<BaseFab<int> >(grids, dmap, 1, 1));
// covered : ghost cells covered by valid cells of this FabArray
// (including periodically shifted valid cells)
// notcovered: ghost cells not covered by valid cells
// (including ghost cells outside periodic boundaries)
// physbnd : boundary cells outside the domain (excluding periodic boundaries)
// interior : interior cells (i.e., valid cells)
int covered = -1;
int interior = 0;
int notcovered = 1; // -> boundary
int physbnd = 2;
masks_m[lev]->BuildMask(geom.Domain(), geom.periodicity(),
covered, notcovered, physbnd, interior);
}
// void buildGeometry_m(const AmrField_t& rho, const AmrField_t& phi,
// const Geometry& geom, int lev, EllipticDomain& bp);
// void solveLevelGeometry_m(AmrField_t& phi,
// const AmrField_t& rho,
// const Geometry& geom, int lev, EllipticDomain& bp);
// void fillMatrixGeometry_m(const Geometry& geom, const AmrField_t& phi, int lev, EllipticDomain& bp);
void setBoundaryValue_m(AmrField_t& phi, const AmrField_t& crse_phi = 0,
IntVect crse_ratio = IntVect::TheZeroVector())
{
// The values of phys_bc & ref_ratio do not matter
// because we are not going to use those parts of MacBndry.
IntVect ref_ratio = IntVect::TheZeroVector();
Array<int> lo_bc(AMREX_SPACEDIM, 0);
Array<int> hi_bc(AMREX_SPACEDIM, 0);
BCRec phys_bc(lo_bc.dataPtr(), hi_bc.dataPtr());
// if (crse_phi == 0 && phi == 0)
// {
// bndry_m.setHomogValues(phys_bc, ref_ratio);
// }
if (crse_phi == 0 && phi != 0)
{
bndry_m->setBndryValues(*phi, 0, 0, phi->nComp(), phys_bc);
}
else if (crse_phi != 0 && phi != 0)
{
BL_ASSERT(crse_ratio != IntVect::TheZeroVector());
const int ncomp = phi->nComp();
const int in_rad = 0;
const int out_rad = 1;
const int extent_rad = 2;
BoxArray crse_boxes(phi->boxArray());
crse_boxes.coarsen(crse_ratio);
BndryRegister crse_br(crse_boxes, phi->DistributionMap(),
in_rad, out_rad, extent_rad, ncomp);
crse_br.copyFrom(*crse_phi, crse_phi->nGrow(), 0, 0, ncomp);
bndry_m->setBndryValues(crse_br, 0, *phi, 0, 0, ncomp, crse_ratio, phys_bc, 3);
// print_m(phi);
}
else
{
amrex::Abort("FMultiGrid::Boundary::build_bndry: How did we get here?");
}
}
void print_m(const AmrField_t& phi) {
for (MFIter mfi(*phi, false); mfi.isValid(); ++mfi) {
const Box& bx = mfi.validbox();
const FArrayBox& fab = (*phi)[mfi];
const int* lo = bx.loVect();
const int* hi = bx.hiVect();
for (int i = lo[0]-1; i <= hi[0]+1; ++i) {
for (int j = lo[1]-1; j <= hi[1]+1; ++j) {
for (int k = lo[2]-1; k <= hi[2]+1; ++k) {
IntVect iv(i, j, k);
std::cout << iv << " " << fab(iv) << std::endl;// std::cin.get();
}
}
}
}
}
void boundaryprint_m() {
bc_m.clear();
for (OrientationIter oitr; oitr; ++oitr) {
const FabSet& fs = bndry_m->bndryValues(oitr());
for (FabSetIter fsi(fs); fsi.isValid(); ++fsi) {
const Box& bx = fsi.validbox();
const FArrayBox& fab = fs[fsi];
const int* lo = bx.loVect();
const int* hi = bx.hiVect();
for (int ii = lo[0]; ii <= hi[0]; ++ii) {
for (int j = lo[1]; j <= hi[1]; ++j) {
for (int k = lo[2]; k <= hi[2]; ++k) {
IntVect iv(ii, j, k);
bc_m.insert(iv);
// std::cout << iv << " " << fab(iv) << std::endl; std::cin.get();
}
}
}
}
}
}
void initGhosts_m() {
for (OrientationIter fi; fi; ++fi)
{
const Orientation face = fi();
FabSet& fs = (*bndry_m.get())[face];
for (FabSetIter fsi(fs); fsi.isValid(); ++fsi)
{
fs[fsi].setVal(0);
}
}
}
void fillDomainBoundary_m(MultiFab& phi, const Geometry& geom) {
Array<BCRec> bc(phi.nComp());
for (int n = 0; n < phi.nComp(); ++n) {
for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
if ( Geometry::isPeriodic(idim) ) {
bc[n].setLo(idim, BCType::int_dir); // interior
bc[n].setHi(idim, BCType::int_dir);
} else {
bc[n].setLo(idim, BCType::foextrap); // first oder extrapolation from last cell in interior
bc[n].setHi(idim, BCType::foextrap);
}
}
}
// fill physical boundary + ghosts at physical boundary
amrex::FillDomainBoundary(phi, geom, bc);
}
void computeFlux_m(const AmrField_t& phi, const Geometry& geom, int lev) {
fluxes_m.resize(lev + 1);
fluxes_m[lev].resize(AMREX_SPACEDIM);
for (int n = 0; n < AMREX_SPACEDIM; ++n) {
BoxArray ba = phi->boxArray();
const DistributionMapping& dmap = phi->DistributionMap();
fluxes_m[lev][n].reset(new MultiFab(ba.surroundingNodes(n), dmap, 1, 1));
fluxes_m[lev][n]->setVal(0.0, 1);
}
const double* dx = geom.CellSize();
for (MFIter mfi(*masks_m[lev], false); mfi.isValid(); ++mfi) {
const Box& bx = mfi.validbox();
const FArrayBox& pfab = (*phi)[mfi];
FArrayBox& xface = (*fluxes_m[lev][0])[mfi];
FArrayBox& yface = (*fluxes_m[lev][1])[mfi];
FArrayBox& zface = (*fluxes_m[lev][2])[mfi];
const int* lo = bx.loVect();
const int* hi = bx.hiVect();
for (int i = lo[0]; i <= hi[0]+1; ++i) {
for (int j = lo[1]; j <= hi[1]+1; ++j) {
for (int k = lo[2]; k <= hi[2]+1; ++k) {
IntVect iv(i, j, k);
// x direction
IntVect left(i-1, j, k);
xface(iv) = (pfab(left) - pfab(iv)) / dx[0];
// y direction
IntVect down(i, j-1, k);
yface(iv) = (pfab(down) - pfab(iv)) / dx[1];
// z direction
IntVect front(i, j, k-1);
zface(iv) = (pfab(front) - pfab(iv)) / dx[2];
}
}
}
}
}
void correctFlux_m(container_t& phi, const Array<Geometry>& geom, int lev) {
const BoxArray& grids = phi[lev+1]->boxArray();
const DistributionMapping& dmap = phi[lev+1]->DistributionMap();
IntVect crse_ratio(2, 2, 2);
std::unique_ptr<FluxRegister> fine = std::unique_ptr<FluxRegister>(new FluxRegister(grids,
dmap,
crse_ratio,
lev+1,
1));
computeFlux_m(phi[lev+1], geom[lev+1], lev+1);
for (int i = 0; i < AMREX_SPACEDIM ; i++) {
Array<MultiFab*> p = amrex::GetArrOfPtrs(fluxes_m[lev+1]);
fine->FineAdd((*p[i]), i, 0, 0, 1, 1.);
}
fine->Reflux(*phi[lev], 1.0, 0, 0, 1, geom[lev]);
}
void fixRHSForSolve(Array<std::unique_ptr<MultiFab> >& rhs,
const Array<std::unique_ptr<FabArray<BaseFab<int> > > >& masks,
const Array<Geometry>& geom, const IntVect& ratio)
{
int num_levels = rhs.size();
for (int lev = 0; lev < num_levels; ++lev) {
MultiFab& fine_rhs = *rhs[lev];
const FabArray<BaseFab<int> >& mask = *masks[lev];
const BoxArray& fine_ba = fine_rhs.boxArray();
const DistributionMapping& fine_dm = fine_rhs.DistributionMap();
MultiFab fine_bndry_data(fine_ba, fine_dm, 1, 1);
zeroOutBoundary(fine_rhs, fine_bndry_data, mask);
}
}
void zeroOutBoundary(MultiFab& input_data,
MultiFab& bndry_data,
const FabArray<BaseFab<int> >& mask)
{
bndry_data.setVal(0.0, 1);
for (MFIter mfi(input_data); mfi.isValid(); ++mfi) {
const Box& bx = mfi.validbox();
zero_out_bndry(bx.loVect(), bx.hiVect(),
input_data[mfi].dataPtr(),
bndry_data[mfi].dataPtr(),
mask[mfi].dataPtr());
}
bndry_data.FillBoundary();
}
private:
AmrField_t fab_set_m;
std::unique_ptr<MacBndry> bndry_m;
Array<std::unique_ptr<FabArray<BaseFab<int> > > > masks_m;
Array<container_t> fluxes_m;
std::set<IntVect> bc_m;
Epetra_MpiComm epetra_comm_m;
Teuchos::RCP<Epetra_Map> map_mp;
int nLevel_m;
// info for a level
IntVect nPoints_m;
// (x, y, z) --> idx
std::map<int, IntVect> indexMap_m;
Teuchos::RCP<Epetra_CrsMatrix> A_mp;
Teuchos::RCP<Epetra_Vector> rho_mp;
Teuchos::RCP<Epetra_Vector> x_mp;
};
#endif
ippl/test/AMR/trilinos/EllipticDomain.cpp deleted 100644 → 0
+ 0
− 667
View file @ 309a1498
// #ifdef HAVE_SAAMG_SOLVER
#include "EllipticDomain.h"
#include <map>
#include <cmath>
#include <iostream>
#include <assert.h>
//FIXME: ORDER HOW TO TRAVERSE NODES IS FIXED, THIS SHOULD BE MORE GENERIC! (PLACES MARKED)
EllipticDomain::EllipticDomain(Vector_t nr, Vector_t hr) {
setNr(nr);
setHr(hr);
}
EllipticDomain::EllipticDomain(double semimajor, double semiminor, Vector_t nr, Vector_t hr, std::string interpl) {
SemiMajor = semimajor;
SemiMinor = semiminor;
setNr(nr);
setHr(hr);
if(interpl == "CONSTANT")
interpolationMethod = CONSTANT;
else if(interpl == "LINEAR")
interpolationMethod = LINEAR;
else if(interpl == "QUADRATIC")
interpolationMethod = QUADRATIC;
}
// EllipticDomain::EllipticDomain(BoundaryGeometry *bgeom, Vector_t nr, Vector_t hr, std::string interpl) {
// SemiMajor = bgeom->getA();
// SemiMinor = bgeom->getB();
// setMinMaxZ(bgeom->getS(), bgeom->getLength());
// setNr(nr);
// setHr(hr);
//
// if(interpl == "CONSTANT")
// interpolationMethod = CONSTANT;
// else if(interpl == "LINEAR")
// interpolationMethod = LINEAR;
// else if(interpl == "QUADRATIC")
// interpolationMethod = QUADRATIC;
// }
//
// EllipticDomain::EllipticDomain(BoundaryGeometry *bgeom, Vector_t nr, std::string interpl) {
// SemiMajor = bgeom->getA();
// SemiMinor = bgeom->getB();
// setMinMaxZ(bgeom->getS(), bgeom->getLength());
// Vector_t hr_m;
// hr_m[0] = (getXRangeMax()-getXRangeMin())/nr[0];
// hr_m[1] = (getYRangeMax()-getYRangeMin())/nr[1];
// hr_m[2] = (getZRangeMax()-getZRangeMin())/nr[2];
// setHr(hr_m);
//
// if(interpl == "CONSTANT")
// interpolationMethod = CONSTANT;
// else if(interpl == "LINEAR")
// interpolationMethod = LINEAR;
// else if(interpl == "QUADRATIC")
// interpolationMethod = QUADRATIC;
// }
EllipticDomain::~EllipticDomain() {
//nothing so far
}
// for this geometry we only have to calculate the intersection on
// one x-y-plane
// for the moment we center the ellipse around the center of the grid
void EllipticDomain::compute(Vector_t hr){
// //there is nothing to be done if the mesh spacings have not changed
// if(hr[0] == getHr()[0] && hr[1] == getHr()[1] && hr[2] == getHr()[2]) {
// hasGeometryChanged_m = false;
// return;
// }
setHr(hr);
hasGeometryChanged_m = true;
//reset number of points inside domain
nxy_m = 0;
// clear previous coordinate maps
IdxMap.clear();
CoordMap.clear();
//clear previous intersection points
IntersectYDir.clear();
IntersectXDir.clear();
// build a index and coordinate map
int idx = 0;
int x, y;
// std::cout << "nr[0] = " << nr[0] << " nr[1] = " << nr[1] << std::endl;
for(x = 0; x < nr[0]; x++) {
for(y = 0; y < nr[1]; y++) {
if(isInside(x, y, 1)) {
IdxMap[toCoordIdx(x, y)] = idx;
CoordMap[idx++] = toCoordIdx(x, y);
// std::cout << idx << " " << toCoordIdx(x, y) << std::endl; std::cin.get();
nxy_m++;
}
}
}
switch(interpolationMethod) {
case CONSTANT:
break;
case LINEAR:
case QUADRATIC:
double smajsq = SemiMajor * SemiMajor;
double sminsq = SemiMinor * SemiMinor;
double yd = 0.0;
double xd = 0.0;
double pos = 0.0;
double mx = (nr[0] - 1) * hr[0] / 2.0;
double my = (nr[1] - 1) * hr[1] / 2.0;
//calculate intersection with the ellipse
for(x = 0; x < nr[0]; x++) {
pos = x * hr[0] - mx;
if (pos <= -SemiMajor || pos >= SemiMajor)
{
IntersectYDir.insert(std::pair<int, double>(x, 0));
IntersectYDir.insert(std::pair<int, double>(x, 0));
}else{
yd = std::abs(sqrt(sminsq - sminsq * pos * pos / smajsq)); // + 0.5*nr[1]*hr[1]);
IntersectYDir.insert(std::pair<int, double>(x, yd));
IntersectYDir.insert(std::pair<int, double>(x, -yd));
}
}
for(y = 0; y < nr[1]; y++) {
pos = y * hr[1] - my;
if (pos <= -SemiMinor || pos >= SemiMinor)
{
IntersectXDir.insert(std::pair<int, double>(y, 0));
IntersectXDir.insert(std::pair<int, double>(y, 0));
}else{
xd = std::abs(sqrt(smajsq - smajsq * pos * pos / sminsq)); // + 0.5*nr[0]*hr[0]);
IntersectXDir.insert(std::pair<int, double>(y, xd));
IntersectXDir.insert(std::pair<int, double>(y, -xd));
}
}
}
}
void EllipticDomain::compute(Vector_t hr, NDIndex<3> localId){
//there is nothing to be done if the mesh spacings have not changed
if(hr[0] == getHr()[0] && hr[1] == getHr()[1] && hr[2] == getHr()[2]) {
hasGeometryChanged_m = false;
return;
}
setHr(hr);
hasGeometryChanged_m = true;
//reset number of points inside domain
nxy_m = 0;
// clear previous coordinate maps
IdxMap.clear();
CoordMap.clear();
//clear previous intersection points
IntersectYDir.clear();
IntersectXDir.clear();
// build a index and coordinate map
int idx = 0;
int x, y;
for(x = localId[0].first(); x<= localId[0].last(); x++) {
for(y = localId[1].first(); y <= localId[1].last(); y++) {
if(isInside(x, y, 1)) {
IdxMap[toCoordIdx(x, y)] = idx;
CoordMap[idx++] = toCoordIdx(x, y);
nxy_m++;
}
}
}
switch(interpolationMethod) {
case CONSTANT:
break;
case LINEAR:
case QUADRATIC:
double smajsq = SemiMajor * SemiMajor;
double sminsq = SemiMinor * SemiMinor;
double yd = 0.0;
double xd = 0.0;
double pos = 0.0;
double mx = (nr[0] - 1) * hr[0] / 2.0;
double my = (nr[1] - 1) * hr[1] / 2.0;
//calculate intersection with the ellipse
for(x = localId[0].first(); x <= localId[0].last(); x++) {
pos = x * hr[0] - mx;
if (pos <= -SemiMajor || pos >= SemiMajor)
{
IntersectYDir.insert(std::pair<int, double>(x, 0));
IntersectYDir.insert(std::pair<int, double>(x, 0));
}else{
yd = std::abs(sqrt(sminsq - sminsq * pos * pos / smajsq)); // + 0.5*nr[1]*hr[1]);
IntersectYDir.insert(std::pair<int, double>(x, yd));
IntersectYDir.insert(std::pair<int, double>(x, -yd));
}
}
for(y = localId[0].first(); y < localId[1].last(); y++) {
pos = y * hr[1] - my;
if (pos <= -SemiMinor || pos >= SemiMinor)
{
IntersectXDir.insert(std::pair<int, double>(y, 0));
IntersectXDir.insert(std::pair<int, double>(y, 0));
}else{
xd = std::abs(sqrt(smajsq - smajsq * pos * pos / sminsq)); // + 0.5*nr[0]*hr[0]);
IntersectXDir.insert(std::pair<int, double>(y, xd));
IntersectXDir.insert(std::pair<int, double>(y, -xd));
}
}
}
}
void EllipticDomain::getBoundaryStencil(int x, int y, int z, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor) {
scaleFactor = 1.0;
// determine which interpolation method we use for points near the boundary
switch(interpolationMethod) {
case CONSTANT:
constantInterpolation(x, y, z, W, E, S, N, F, B, C, scaleFactor);
break;
case LINEAR:
linearInterpolation(x, y, z, W, E, S, N, F, B, C, scaleFactor);
break;
case QUADRATIC:
quadraticInterpolation(x, y, z, W, E, S, N, F, B, C, scaleFactor);
break;
}
// stencil center value has to be positive!
assert(C > 0);
}
void EllipticDomain::getBoundaryStencil(int idx, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor) {
int x = 0, y = 0, z = 0;
getCoord(idx, x, y, z);
getBoundaryStencil(x, y, z, W, E, S, N, F, B, C, scaleFactor);
}
void EllipticDomain::getNeighbours(int idx, int &W, int &E, int &S, int &N, int &F, int &B) {
int x = 0, y = 0, z = 0;
getCoord(idx, x, y, z);
getNeighbours(x, y, z, W, E, S, N, F, B);
}
void EllipticDomain::getNeighbours(int x, int y, int z, int &W, int &E, int &S, int &N, int &F, int &B) {
if(x > 0)
W = getIdx(x - 1, y, z);
else
W = -1;
if(x < nr[0] - 1)
E = getIdx(x + 1, y, z);
else
E = -1;
if(y < nr[1] - 1)
N = getIdx(x, y + 1, z);
else
N = -1;
if(y > 0)
S = getIdx(x, y - 1, z);
else
S = -1;
if(z > 0)
F = getIdx(x, y, z - 1);
else
F = -1;
if(z < nr[2] - 1)
B = getIdx(x, y, z + 1);
else
B = -1;
}
void EllipticDomain::constantInterpolation(int x, int y, int z, double &WV, double &EV, double &SV, double &NV, double &FV, double &BV, double &CV, double &scaleFactor) {
scaleFactor = 1.0;
WV = -1/(hr[0]*hr[0]);
EV = -1/(hr[0]*hr[0]);
NV = -1/(hr[1]*hr[1]);
SV = -1/(hr[1]*hr[1]);
FV = -1/(hr[2]*hr[2]);
BV = -1/(hr[2]*hr[2]);
CV = 2/(hr[0]*hr[0]) + 2/(hr[1]*hr[1]) + 2/(hr[2]*hr[2]);
// we are a right boundary point
if(!isInside(x + 1, y, z))
EV= 0.0;
// we are a left boundary point
if(!isInside(x - 1, y, z))
WV = 0.0;
// we are a upper boundary point
if(!isInside(x, y + 1, z))
NV = 0.0;
// we are a lower boundary point
if(!isInside(x, y - 1, z))
SV = 0.0;
if(z == 1 || z == nr[2] - 2) {
// case where we are on the Robin BC in Z-direction
// where we distinguish two cases
// IFF: this values should not matter because they
// never make it into the discretization matrix
if(z == 1)
FV = 0.0;
else
BV = 0.0;
// add contribution of Robin discretization to center point
// d the distance between the center of the bunch and the boundary
//double cx = (x-(nr[0]-1)/2)*hr[0];
//double cy = (y-(nr[1]-1)/2)*hr[1];
//double cz = hr[2]*(nr[2]-1);
//double d = sqrt(cx*cx+cy*cy+cz*cz);
double d = hr[2] * (nr[2] - 1) / 2;
CV += 2 / (d * hr[2]);
//C += 2/((hr[2]*(nr[2]-1)/2.0) * hr[2]);
// scale all stencil-points in z-plane with 0.5 (Robin discretization)
WV /= 2.0;
EV /= 2.0;
NV /= 2.0;
SV /= 2.0;
CV /= 2.0;
scaleFactor *= 0.5;
}
}
void EllipticDomain::linearInterpolation(int x, int y, int z, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor) {
scaleFactor = 1.0;
double cx = x * hr[0] - (nr[0] - 1) * hr[0] / 2.0;
double cy = y * hr[1] - (nr[1] - 1) * hr[1] / 2.0;
double dx = 0.0;
std::multimap<int, double>::iterator it = IntersectXDir.find(y);
if(cx < 0)
it++;
dx = it->second;
double dy = 0.0;
it = IntersectYDir.find(x);
if(cy < 0)
it++;
dy = it->second;
double dw = hr[0];
double de = hr[0];
double dn = hr[1];
double ds = hr[1];
C = 0.0;
//we are a right boundary point
if(!isInside(x + 1, y, z)) {
C += 1 / ((dx - cx) * de);
E = 0.0;
} else {
C += 1 / (de * de);
E = -1 / (de * de);
}
//we are a left boundary point
if(!isInside(x - 1, y, z)) {
C += 1 / ((std::abs(dx) - std::abs(cx)) * dw);
W = 0.0;
} else {
C += 1 / (dw * dw);
W = -1 / (dw * dw);
}
//we are a upper boundary point
if(!isInside(x, y + 1, z)) {
C += 1 / ((dy - cy) * dn);
N = 0.0;
} else {
C += 1 / (dn * dn);
N = -1 / (dn * dn);
}
//we are a lower boundary point
if(!isInside(x, y - 1, z)) {
C += 1 / ((std::abs(dy) - std::abs(cy)) * ds);
S = 0.0;
} else {
C += 1 / (ds * ds);
S = -1 / (ds * ds);
}
F = -1 / (hr[2] * hr[2]);
B = -1 / (hr[2] * hr[2]);
C += 2 / (hr[2] * hr[2]);
// handle boundary condition in z direction
/*
if(z == 0 || z == nr[2] - 1) {
// case where we are on the NEUMAN BC in Z-direction
// where we distinguish two cases
if(z == 0)
F = 0.0;
else
B = 0.0;
//hr[2]*(nr2[2]-1)/2 = radius
double d = hr[2] * (nr[2] - 1) / 2;
C += 2 / (d * hr[2]);
W /= 2.0;
E /= 2.0;
N /= 2.0;
S /= 2.0;
C /= 2.0;
scaleFactor *= 0.5;
}
*/
}
void EllipticDomain::quadraticInterpolation(int x, int y, int z, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor) {
double cx = (x - (nr[0] - 1) / 2.0) * hr[0];
double cy = (y - (nr[1] - 1) / 2.0) * hr[1];
// since every vector for elliptic domains has ALWAYS size 2 we
// can catch all cases manually
double dx = 0.0;
std::multimap<int, double>::iterator it = IntersectXDir.find(y);
if(cx < 0)
it++;
dx = it->second;
double dy = 0.0;
it = IntersectYDir.find(x);
if(cy < 0)
it++;
dy = it->second;
double dw = hr[0];
double de = hr[0];
double dn = hr[1];
double ds = hr[1];
W = 1.0;
E = 1.0;
N = 1.0;
S = 1.0;
F = 1.0;
B = 1.0;
C = 0.0;
//TODO: = cx+hr[0] > dx && cx > 0
//if((x-nr[0]/2.0+1)*hr[0] > dx && cx > 0) {
////we are a right boundary point
////if(!isInside(x+1,y,z)) {
//de = dx-cx;
//}
//if((x-nr[0]/2.0-1)*hr[0] < dx && cx < 0) {
////we are a left boundary point
////if(!isInside(x-1,y,z)) {
//dw = std::abs(dx)-std::abs(cx);
//}
//if((y-nr[1]/2.0+1)*hr[1] > dy && cy > 0) {
////we are a upper boundary point
////if(!isInside(x,y+1,z)) {
//dn = dy-cy;
//}
//if((y-nr[1]/2.0-1)*hr[1] < dy && cy < 0) {
////we are a lower boundary point
////if(!isInside(x,y-1,z)) {
//ds = std::abs(dy)-std::abs(cy);
//}
//TODO: = cx+hr[0] > dx && cx > 0
//if((x-nr[0]/2.0+1)*hr[0] > dx && cx > 0) {
//we are a right boundary point
if(!isInside(x + 1, y, z)) {
de = dx - cx;
E = 0.0;
}
//if((x-nr[0]/2.0-1)*hr[0] < dx && cx < 0) {
//we are a left boundary point
if(!isInside(x - 1, y, z)) {
dw = std::abs(dx) - std::abs(cx);
W = 0.0;
}
//if((y-nr[1]/2.0+1)*hr[1] > dy && cy > 0) {
//we are a upper boundary point
if(!isInside(x, y + 1, z)) {
dn = dy - cy;
N = 0.0;
}
//if((y-nr[1]/2.0-1)*hr[1] < dy && cy < 0) {
//we are a lower boundary point
if(!isInside(x, y - 1, z)) {
ds = std::abs(dy) - std::abs(cy);
S = 0.0;
}
//2/dw*(dw_de)
W *= -1.0 / (dw * (dw + de));
E *= -1.0 / (de * (dw + de));
N *= -1.0 / (dn * (dn + ds));
S *= -1.0 / (ds * (dn + ds));
F = -1 / (hr[2] * (hr[2] + hr[2]));
B = -1 / (hr[2] * (hr[2] + hr[2]));
//TODO: problem when de,dw,dn,ds == 0
//is NOT a regular BOUND PT
C += 1 / de * 1 / dw;
C += 1 / dn * 1 / ds;
C += 1 / hr[2] * 1 / hr[2];
scaleFactor = 0.5;
//for regular gridpoints no problem with symmetry, just boundary
//z direction is right
//implement isLastInside(dir)
//we have LOCAL x,y coordinates!
/*
if(dw != 0 && !wIsB)
W = -1/dw * (dn+ds) * 2*hr[2];
else
W = 0;
if(de != 0 && !eIsB)
E = -1/de * (dn+ds) * 2*hr[2];
else
E = 0;
if(dn != 0 && !nIsB)
N = -1/dn * (dw+de) * 2*hr[2];
else
N = 0;
if(ds != 0 && !sIsB)
S = -1/ds * (dw+de) * 2*hr[2];
else
S = 0;
F = -(dw+de)*(dn+ds)/hr[2];
B = -(dw+de)*(dn+ds)/hr[2];
*/
//if(dw != 0)
//W = -2*hr[2]*(dn+ds)/dw;
//else
//W = 0;
//if(de != 0)
//E = -2*hr[2]*(dn+ds)/de;
//else
//E = 0;
//if(dn != 0)
//N = -2*hr[2]*(dw+de)/dn;
//else
//N = 0;
//if(ds != 0)
//S = -2*hr[2]*(dw+de)/ds;
//else
//S = 0;
//F = -(dw+de)*(dn+ds)/hr[2];
//B = -(dw+de)*(dn+ds)/hr[2];
//// RHS scaleFactor for current 3D index
//// Factor 0.5 results from the SW/quadratic extrapolation
//scaleFactor = 0.5*(dw+de)*(dn+ds)*(2*hr[2]);
// catch the case where a point lies on the boundary
//FIXME: do this more elegant!
//double m1 = dw*de;
//double m2 = dn*ds;
//if(de == 0)
//m1 = dw;
//if(dw == 0)
//m1 = de;
//if(dn == 0)
//m2 = ds;
//if(ds == 0)
//m2 = dn;
////XXX: dn+ds || dw+de can be 0
////C = 2*(dn+ds)*(dw+de)/hr[2];
//C = 2/hr[2];
//if(dw != 0 || de != 0)
//C *= (dw+de);
//if(dn != 0 || ds != 0)
//C *= (dn+ds);
//if(dw != 0 || de != 0)
//C += (2*hr[2])*(dn+ds)*(dw+de)/m1;
//if(dn != 0 || ds != 0)
//C += (2*hr[2])*(dw+de)*(dn+ds)/m2;
//handle Neumann case
//if(z == 0 || z == nr[2]-1) {
//if(z == 0)
//F = 0.0;
//else
//B = 0.0;
////neumann stuff
//W = W/2.0;
//E = E/2.0;
//N = N/2.0;
//S = S/2.0;
//C /= 2.0;
//scaleFactor /= 2.0;
//}
// handle boundary condition in z direction
if(z == 0 || z == nr[2] - 1) {
// case where we are on the NEUMAN BC in Z-direction
// where we distinguish two cases
if(z == 0)
F = 0.0;
else
B = 0.0;
//C += 2/((hr[2]*(nr[2]-1)/2.0) * hr[2]);
//hr[2]*(nr2[2]-1)/2 = radius
double d = hr[2] * (nr[2] - 1) / 2;
C += 2 / (d * hr[2]);
W /= 2.0;
E /= 2.0;
N /= 2.0;
S /= 2.0;
C /= 2.0;
scaleFactor /= 2.0;
}
}
// #endif //#ifdef HAVE_SAAMG_SOLVER
ippl/test/AMR/trilinos/EllipticDomain.h deleted 100644 → 0
+ 0
− 119
View file @ 309a1498
#ifndef ELLIPTICAL_DOMAIN_H
#define ELLIPTICAL_DOMAIN_H
// #ifdef HAVE_SAAMG_SOLVER
#include <vector>
#include <map>
#include <string>
#include <cmath>
#include "IrregularDomain.h"
// #include "Structure/BoundaryGeometry.h"
class EllipticDomain : public IrregularDomain {
public:
EllipticDomain(Vector_t nr, Vector_t hr);
EllipticDomain(double semimajor, double semiminor, Vector_t nr, Vector_t hr, std::string interpl);
// EllipticDomain(BoundaryGeometry *bgeom, Vector_t nr, Vector_t hr, std::string interpl);
// EllipticDomain(BoundaryGeometry *bgeom, Vector_t nr, std::string interpl);
~EllipticDomain();
/// returns discretization at (x,y,z)
void getBoundaryStencil(int x, int y, int z, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor);
/// returns discretization at 3D index
void getBoundaryStencil(int idx, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor);
/// returns index of neighbours at (x,y,z)
void getNeighbours(int x, int y, int z, int &W, int &E, int &S, int &N, int &F, int &B);
/// returns index of neighbours at 3D index
void getNeighbours(int idx, int &W, int &E, int &S, int &N, int &F, int &B);
/// returns type of boundary condition
std::string getType() {return "Elliptic";}
/// queries if a given (x,y,z) coordinate lies inside the domain
inline bool isInside(int x, int y, int z) {
double xx = (x - (nr[0] - 1) / 2.0) * hr[0];
double yy = (y - (nr[1] - 1) / 2.0) * hr[1];
return ((xx * xx / (SemiMajor * SemiMajor) + yy * yy / (SemiMinor * SemiMinor) < 1) && z != 0 && z != nr[2] - 1);
}
int getNumXY(int z) { return nxy_m; }
/// set semi-minor
void setSemiMinor(double sm) {SemiMinor = sm;}
/// set semi-major
void setSemiMajor(double sm) {SemiMajor = sm;}
/// calculates intersection
void compute(Vector_t hr);
void compute(Vector_t hr, NDIndex<3> localId);
double getXRangeMin() { return -SemiMajor; }
double getXRangeMax() { return SemiMajor; }
double getYRangeMin() { return -SemiMinor; }
double getYRangeMax() { return SemiMinor; }
double getZRangeMin() { return zMin_m; }
double getZRangeMax() { return zMax_m; }
//TODO: ?
int getStartIdx() {return 0;}
bool hasGeometryChanged() { return hasGeometryChanged_m; }
private:
/// Map from a single coordinate (x or y) to a list of intersection values with
/// boundary.
typedef std::multimap<int, double> EllipticPointList;
/// all intersection points with grid lines in X direction
EllipticPointList IntersectXDir;
/// all intersection points with grid lines in Y direction
EllipticPointList IntersectYDir;
/// mapping (x,y,z) -> idx
std::map<int, int> IdxMap;
/// mapping idx -> (x,y,z)
std::map<int, int> CoordMap;
/// semi-major of the ellipse
double SemiMajor;
/// semi-minor of the ellipse
double SemiMinor;
/// number of nodes in the xy plane (for this case: independent of the z coordinate)
int nxy_m;
/// interpolation type
int interpolationMethod;
/// flag indicating if geometry has changed for the current time-step
bool hasGeometryChanged_m;
/// conversion from (x,y) to index in xy plane
inline int toCoordIdx(int x, int y) { return y * nr[0] + x; }
/// conversion from (x,y,z) to index on the 3D grid
inline int getIdx(int x, int y, int z) {
if(isInside(x, y, z) && x >= 0 && y >= 0 && z >= 0)
return IdxMap[toCoordIdx(x, y)] + (z - 1) * nxy_m;
else
return -1;
}
/// conversion from a 3D index to (x,y,z)
inline void getCoord(int idx, int &x, int &y, int &z) {
int ixy = idx % nxy_m;
int xy = CoordMap[ixy];
int inr = (int)nr[0];
x = xy % inr;
y = (xy - x) / nr[0];
z = (idx - ixy) / nxy_m + 1;
}
/// different interpolation methods for boundary points
void constantInterpolation(int x, int y, int z, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor);
void linearInterpolation(int x, int y, int z, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor);
void quadraticInterpolation(int x, int y, int z, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor);
};
// #endif //#ifdef HAVE_SAAMG_SOLVER
#endif //#ifdef ELLIPTICAL_DOMAIN_H
ippl/test/AMR/trilinos/IrregularDomain.h deleted 100644 → 0
+ 0
− 129
View file @ 309a1498
#ifndef IRREGULAR_DOMAIN_H
#define IRREGULAR_DOMAIN_H
// #ifdef HAVE_SAAMG_SOLVER
#include "Ippl.h"
#include <vector>
#include <string>
// #include "Algorithms/PBunchDefs.h"
// #include "Algorithms/Quaternion.h"
typedef Vektor<double, AMREX_SPACEDIM> Vector_t;
/// enumeration corresponding to different interpolation methods at the boundary
enum {
CONSTANT,
LINEAR,
QUADRATIC
};
/// this class defines a common abstract interface for different types of boundaries
class IrregularDomain {
public:
/** method to compute the intersection points with the boundary geometry (stored in some appropriate data structure)
* \param hr updated mesh spacings
*/
virtual void compute(Vector_t hr) = 0;
virtual void compute(Vector_t hr, NDIndex<3> localId) = 0;
/** method to get the number of gridpoints in a given z plane
* \param z coordinate of the z plane
* \return int number of grid nodes in the given z plane
*/
virtual int getNumXY(int z) = 0;
/// method to calculate the stencil at a boundary points
/// \param x index of the current element in the matrix
/// \param y index of the current element in the matrix
/// \param z index of the current element in the matrix
/// \param W stencil value of the element in the west of idx: (x-1)
/// \param E stencil value of the element in the east of idx: (x+1)
/// \param S stencil value of the element in the south of idx: (y-1)
/// \param N stencil value of the element in the north of idx: (y+1)
/// \param F stencil value of the element in front of idx: (z-1)
/// \param B stencil value of the element in the back of idx: (z+1)
/// \param C stencil value of the element in the center
virtual void getBoundaryStencil(int x, int y, int z, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor) = 0;
/// method to calculate the stencil at a boundary points
/// \param idx index of the current element in the matrix
/// \param W stencil value of the element in the west of idx: (x-1)
/// \param E stencil value of the element in the east of idx: (x+1)
/// \param S stencil value of the element in the south of idx: (y-1)
/// \param N stencil value of the element in the north of idx: (y+1)
/// \param F stencil value of the element in front of idx: (z-1)
/// \param B stencil value of the element in the back of idx: (z+1)
/// \param C stencil value of the element in the center
virtual void getBoundaryStencil(int idx, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor) = 0;
/// method to calculate the neighbours in the matrix of the current index (x,y,z)
/// \param x index of the current element in the matrix
/// \param y index of the current element in the matrix
/// \param z index of the current element in the matrix
/// \param W stencil index of the element in the west of idx: (x-1)
/// \param E stencil index of the element in the east of idx: (x+1)
/// \param S stencil index of the element in the south of idx: (y-1)
/// \param N stencil index of the element in the north of idx: (y+1)
/// \param F stencil index of the element in front of idx: (z-1)
/// \param B stencil index of the element in the back of idx: (z+1)
virtual void getNeighbours(int x, int y, int z, int &W, int &E, int &S, int &N, int &F, int &B) = 0;
virtual void getNeighbours(int idx, int &W, int &E, int &S, int &N, int &F, int &B) = 0;
/// method that identifies a specialized boundary geometry
/// \return std::string containing a description of the boundary geometry used
virtual std::string getType() = 0;
/// method that checks if a given point lies inside the boundary
/// \param x index of the current element in the matrix
/// \param y index of the current element in the matrix
/// \param z index of the current element in the matrix
/// \return boolean indicating if the point lies inside the boundary
virtual bool isInside(int x, int y, int z) = 0;
Vector_t getNr() { return nr; }
Vector_t getHr() { return hr; }
void setNr(Vector_t nri) { nr = nri; }
void setHr(Vector_t hri) { hr = hri; }
void setMinMaxZ(double minz, double maxz) { zMin_m=minz; zMax_m=maxz; }
double getMinZ() { return zMin_m; }
double getMaxZ() { return zMax_m; }
void setGlobalMeanR(Vector_t rmean) { rMean_m = rmean;}
Vector_t getGlobalMeanR() { return rMean_m; }
// void setGlobalToLocalQuaternion(Quaternion_t globalToLocalQuaternion){
// globalToLocalQuaternion_m = globalToLocalQuaternion;}
// Quaternion_t getGlobalToLocalQuaternion() { return globalToLocalQuaternion_m;}
virtual double getXRangeMin() = 0;
virtual double getXRangeMax() = 0;
virtual double getYRangeMin() = 0;
virtual double getYRangeMax() = 0;
virtual double getZRangeMin() = 0;
virtual double getZRangeMax() = 0;
virtual int getIdx(int x, int y, int z) = 0;
virtual bool hasGeometryChanged() = 0;
protected:
// a irregular domain is always defined on a grid
/// number of mesh points in each direction
Vector_t nr;
/// mesh-spacings in each direction
Vector_t hr;
/// min/max of bunch in floor coordinates
double zMin_m;
double zMax_m;
/// mean position of bunch (m)
Vector_t rMean_m;
// Quaternion_t globalToLocalQuaternion_m;
};
// #endif //#ifdef HAVE_SAAMG_SOLVER
#endif //#ifndef IRREGULAR_DOMAIN_H
\ No newline at end of file
ippl/test/AMR/trilinos/testTrilinosGeometry.cpp deleted 100644 → 0
+ 0
− 607
View file @ 309a1498
This diff is collapsed.
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment