Commit 56863f44 authored by frey_m's avatar frey_m
Browse files

Adding AMR solver to OPAL

parent f371a233
......@@ -119,6 +119,7 @@ IF (USE_H5HUT2)
ENDIF (USE_H5HUT2)
IF (ENABLE_AMR_SOLVER)
MESSAGE (STATUS "Enable AMR_SOLVER " ${ENABLE_AMR_SOLVER})
ENABLE_LANGUAGE (Fortran)
SET ( CCSE_DIR $ENV{BOXLIB_PREFIX} )
......@@ -169,9 +170,8 @@ STRING(REGEX MATCH "[^0-9]*" HOSTNAME_BASE "${HOSTNAME}")
# Trilinos_TPL_LIBRARIES, instead of cray-tpsl/16.07.1 it uses
# cray-tpsl/16.06.1 libraries -> Linker error.
# We can avoid this issue by not going into the if-statement
if (NOT ${HOSTNAME_BASE} MATCHES "edison" AND NOT ${HOSTNAME_BASE} MATCHES "cori" AND (ENABLE_SAAMG_SOLVER OR ENABLE_AMR_SOLVER) )
if (NOT ${HOSTNAME_BASE} MATCHES "edison" AND NOT ${HOSTNAME_BASE} MATCHES "cori" AND ENABLE_SAAMG_SOLVER )
MESSAGE (STATUS "Enable SAAMG_SOLVER " ${ENABLE_SAAMG_SOLVER})
MESSAGE (STATUS "Enable AMR_SOLVER " ${ENABLE_AMR_SOLVER})
find_package (Trilinos REQUIRED HINTS $ENV{TRILINOS_PREFIX} $ENV{TRILINOS_DIR} $ENV{TRILINOS})
message (STATUS "Found Trilinos: ${Trilinos_DIR}")
......@@ -194,6 +194,7 @@ if (NOT ${HOSTNAME_BASE} MATCHES "edison" AND NOT ${HOSTNAME_BASE} MATCHES "cori
endif ()
endif ()
IF (DBG_SCALARFIELD)
MESSAGE (STATUS "\nWrite scalar rho_m field is enabled ")
SET (CMAKE_CXX_FLAGS "-DDBG_SCALARFIELD ${CMAKE_CXX_FLAGS}")
......
......@@ -143,7 +143,7 @@ IF (ENABLE_AMR_SOLVER)
MESSAGE (STATUS " ENABLE_PROFILING = ${ENABLE_PROFILING} (INT: 0, 1)")
ENDIF (ENABLE_AMR_SOLVER)
if (ENABLE_SAAMG_SOLVER OR ENABLE_AMR_SOLVER)
if (ENABLE_SAAMG_SOLVER)
find_package (Trilinos REQUIRED HINTS $ENV{TRILINOS_PREFIX} $ENV{TRILINOS_DIR} $ENV{TRILINOS})
message (STATUS "Found Trilinos: ${Trilinos_DIR}")
......
......@@ -15,6 +15,11 @@
#include "Field/Field.h"
#include "FFT/FFT.h"
#ifdef HAVE_AMR_SOLVER
#include <MultiFab.H>
#include <PArray.H>
#endif
typedef IntCIC IntrplCIC_t;
typedef IntNGP IntrplNGP_t;
typedef IntSUDS IntrplSUDS_t;
......@@ -49,6 +54,13 @@ typedef FFT<RCTransform, 3, double> FFT_t;
typedef FFT<SineTransform, 3, double> SINE_t;
typedef FFT<CCTransform, 3, double> FFTC_t;
#ifdef HAVE_AMR_SOLVER
//TODO make it AMR object dependent (i.e. typedef in AMR object class
typedef MultiFab AmrField_t;
typedef PArray<AmrField_t> AmrFieldContainer_t;
#endif
namespace ParticleType {
enum { REGULAR,
FIELDEMISSION,
......
#ifndef AMR_POISSON_SOLVER_H_
#define AMR_POISSON_SOLVER_H_
#include "PoissonSolver.h"
#include <memory>
template<class AmrObject>
class AmrPoissonSolver : public PoissonSolver {
public:
/*!
* @param amrobject_p holds information about grids and domain
*/
AmrPoissonSolver(AmrObject* amrobject_p) : amrobject_mp(amrobject_p) {}
virtual ~AmrPoissonSolver();
void test(PartBunch &bunch) { };
private:
std::unique_ptr<AmrObject> amrobject_mp;
};
#endif
\ No newline at end of file
set (_SRCS
FMGPoissonSolver.cpp
)
include_directories (
${CMAKE_CURRENT_SOURCE_DIR}
)
add_opal_sources(${_SRCS})
set (HDRS
FMGPoissonSolver.h
)
install (FILES ${HDRS} DESTINATION "${CMAKE_INSTALL_PREFIX}/include/Solvers/BoxLibSolvers")
\ No newline at end of file
#include "FMGPoissonSolver.h"
FMGPoissonSolver::FMGPoissonSolver(AmrBoxLib* amrobject_p)
: AmrPoissonSolver<AmrBoxlib>(amrobject_p),
tol_m(1.e-10),
abstol_m(1.e-14),
phi_m(PArrayManage)
{
// Dirichlet boundary conditions are default
for (int d = 0; d < BL_SPACEDIM; ++d) {
bc_m[2 * d] = MGT_BC_DIR;
bc_m[2 * d + 1] = MGT_BC_DIR;
}
}
void FMGPoissonSolver::solve(AmrFieldContainer_t &rho,
AmrFieldContainer_t &efield,
unsigned short baseLevel,
unsigned short finestLevel)
{
const GeomContainer_t& geom = amrobject_p->Geom();
if (Geometry::isAllPeriodic()) {
for (int d = 0; d < BL_SPACEDIM; ++d) {
bc_m[2 * d] = MGT_BC_PER;
bc_m[2 * d + 1] = MGT_BC_PER;
}
} else if ( Geometry::isAnyPeriodic() ) {
for (int d = 0; d < BL_SPACEDIM; ++d) {
if ( Geometry::isPeriodic(d) ) {
bc_m[2 * d] = MGT_BC_PER;
bc_m[2 * d + 1] = MGT_BC_PER;
}
}
}
Array<AmrFieldContainer_t> grad_phi_edge(rho.size(), PArrayManage);
for (int lev = baseLevel; lev <= finestLevel ; ++lev) {
grad_phi_edge[lev].resize(BL_SPACEDIM);
for (int n = 0; n < BL_SPACEDIM; ++n) {
BoxArray ba = rho[lev].boxArray();
grad_phi_edge[lev].set(n, new MultiFab(ba.surroundingNodes(n), 1, 1));
}
}
// initialize potential and electric field grids on each level
initGrids_m(efield);
this->solveWithF90_m(rho, grad_phi_edge, geom,
baseLevel, finestLevel, tol_m, abstol_m);
for (int lev = baseLevel; lev <= finestLevel; ++lev) {
BoxLib::average_face_to_cellcenter(efield[lev],
grad_phi_edge[lev],
geom[lev]);
efield[lev].FillBoundary(0,BL_SPACEDIM,geom[lev].periodicity());
efield[lev].mult(-1.0, 0, 3);
}
}
void FMGPoissonSolver::getXRangeMin(unsigned short level) {
return amrobject_mp->Geom(level).ProbLo(0);
}
void FMGPoissonSolver::getXRangeMax(unsigned short level) {
return amrobject_mp->Geom(level).ProbHi(0);
}
void FMGPoissonSolver::getYRangeMin(unsigned short level) {
return amrobject_mp->Geom(level).ProbLo(1);
}
void FMGPoissonSolver::getYRangeMax(unsigned short level) {
return amrobject_mp->Geom(level).ProbHi(1);
}
void FMGPoissonSolver::getYRangeMin(unsigned short level) {
return amrobject_mp->Geom(level).ProbLo(2);
}
void FMGPoissonSolver::getZRangeMax(unsigned short level) {
return amrobject_mp->Geom(level).ProbHi(2);
}
Inform &FMGPoissonSolver::print(Inform &os) const {
os << "* ************* F M G P o i s s o n S o l v e r ************************************ " << endl
<< "* tolerance " << tol_m << '\n'
<< "* absolute tolerance " << abstol_m << '\n' << endl
<< "* ******************************************************************** " << endl;
return os;
}
void FMGPoissonSolver::solveWithF90_m(AmrFieldContainer_t& rho,
Array<AmrFieldContainer_t>& grad_phi_edge,
const GeomContainer_t& geom,
int baseLevel,
int finestLevel)
{
int nlevs = finestLevel - baseLevel + 1;
PArray<Geometry> geom_p(nlevs, PArrayManage);
AmrFieldContainer_t rho_p(nlevs, PArrayManage);
AmrFieldContainer_t phi_p(nlevs, PArrayManage);
for (int ilev = 0; ilev < nlevs; ++ilev) {
geom_p.set(ilev, &geom[ilev + baseLevel]);
rho_p.set(ilev, &rho[ilev + baseLevel]);
phi_p.set(ilev, &phi_m[ilev + baseLevel]);
}
// Refinement ratio is hardwired to 2 here.
IntVect crse_ratio = (baseLevel == 0) ?
IntVect::TheZeroVector() : IntVect::TheUnitVector() * 2;
FMultiGrid fmg(geom_p, baseLevel, crse_ratio);
if (baseLevel == 0)
fmg.set_bc(bc_m, phi_m[baseLevel]);
else
fmg.set_bc(bc_m, phi_m[baseLevel-1], phi_m[baseLevel]);
/* (alpha * a - beta * (del dot b grad)) phi = rho
* (b = 1)
*
* The function call set_const_gravity_coeffs() sets alpha = 0.0
* and beta = -1 (in MGT_Solver::set_const_gravity_coeffs)
*
* --> (del dot grad) phi = rho
*/
fmg.set_const_gravity_coeffs();
fmg.set_maxorder(3);
int always_use_bnorm = 0;
int need_grad_phi = 1;
fmg.solve(phi_p, rho_p, tol_m, abstol_m, always_use_bnorm, need_grad_phi);
for (int ilev = 0; ilev < nlevs; ++ilev) {
int amr_level = ilev + baseLevel;
fmg.get_fluxes(grad_phi_edge[amr_level], ilev);
}
}
void FMGPoissonSolver::initGrids_m(AmrFieldContainer_t& efield) {
for (std::size_t lev = 0; lev < efield.size(); ++lev) {
phi_m.clear(lev);
efield.clear(lev);
const BoxArray& ba = amrobject_mp->boxArray()[lev];
// #components #ghost cells
phi_m.set(lev, new AmrField_t(ba, 1, 1) );
efield.set(lev, new AmrField_t(ba, 3, 1) );
phi_m[lev].setVal(0.0);
efield[lev].setVal(0.0);
}
}
#ifndef FMG_POISSON_SOLVER_H_
#define FMG_POISSON_SOLVER_H_
#include <BoxLib.H>
#include <MultiFabUtil.H>
#include <BLFort.H>
#include <MacBndry.H>
#include <MGT_Solver.H>
#include <mg_cpp_f.h>
#include <stencil_types.H>
#include <FMultiGrid.H>
class AmrBoxlib; //TODO remove template and do AmrPoissonSolver<AmrBoxlib>
class FMGPoissonSolver : public AmrPoissonSolver<AmrBoxlib> {
private:
typedef Array<Geometry> GeomContainer_t;
public:
FMGPoissonSolver(AmrBoxlib* amrobject_p);
Inform &print(Inform &os) const;
private:
void solveWithF90_m(AmrFieldContainer_t& rho,
Array< AmrFieldContainer_t >& grad_phi_edge,
const GeomContainer_t& geom,
int baseLevel,
int finestLevel);
/*! Initialize the potential and electric field to zero
* on each level and grid.
* @param efield to be intialized.
*/
void initGrids_m(AmrFieldContainer_t& efield);
private:
int bc_m[2*BL_SPACEDIM]; ///< Boundary conditions
double tol_m;
double abstol_m;
AmrFieldContainer_t phi_m;
};
inline Inform &operator<<(Inform &os, const FMGPoissonSolver &fs) {
return fs.print(os);
}
#endif
\ No newline at end of file
......@@ -30,4 +30,10 @@ set (HDRS
TaperDomain.h
)
if ( ENABLE_AMR_SOLVER )
list(APPEND HDRS AmrPoissonSolver.h)
add_subdirectory(BoxLibSolvers)
endif (ENABLE_AMR_SOLVER )
install (FILES ${HDRS} DESTINATION "${CMAKE_INSTALL_PREFIX}/include/Solvers")
\ No newline at end of file
......@@ -54,12 +54,12 @@ public:
/// compute the shifted integrated Green function as described in <A HREF="http://prst-ab.aps.org/abstract/PRSTAB/v9/i4/e044204">Three-dimensional quasistatic model for high brightness beam dynamics simulation</A> by Qiang et al.
void shiftedIntGreensFunction(double zshift);
double getXRangeMin() {return -a_m;}
double getXRangeMax() {return a_m;}
double getYRangeMin() {return -a_m;}
double getYRangeMax() {return a_m;}
double getZRangeMin() {return -a_m; }
double getZRangeMax() {return a_m; }
double getXRangeMin(unsigned short level) {return -a_m;}
double getXRangeMax(unsigned short level) {return a_m;}
double getYRangeMin(unsigned short level) {return -a_m;}
double getYRangeMax(unsigned short level) {return a_m;}
double getZRangeMin(unsigned short level) {return -a_m; }
double getZRangeMax(unsigned short level) {return a_m; }
void test(PartBunch &bunch) { }
......
......@@ -77,12 +77,12 @@ public:
/// compute the shifted integrated Green function as described in <A HREF="http://prst-ab.aps.org/abstract/PRSTAB/v9/i4/e044204">Three-dimensional quasistatic model for high brightness beam dynamics simulation</A> by Qiang et al.
void shiftedIntGreensFunction(double zshift);
double getXRangeMin() {return 1.0;}
double getXRangeMax() {return 1.0;}
double getYRangeMin() {return 1.0;}
double getYRangeMax() {return 1.0;}
double getZRangeMin() {return 1.0;}
double getZRangeMax() {return 1.0;}
double getXRangeMin(unsigned short level) {return 1.0;}
double getXRangeMax(unsigned short level) {return 1.0;}
double getYRangeMin(unsigned short level) {return 1.0;}
double getYRangeMax(unsigned short level) {return 1.0;}
double getZRangeMin(unsigned short level) {return 1.0;}
double getZRangeMax(unsigned short level) {return 1.0;}
void test(PartBunch &bunch) { }
Inform &print(Inform &os) const;
......
......@@ -117,12 +117,12 @@ public:
/// force Solver to recompute Epetra_Map
void recomputeMap() { hasParallelDecompositionChanged_m = true; }
double getXRangeMin() { return bp->getXRangeMin(); }
double getXRangeMax() { return bp->getXRangeMax(); }
double getYRangeMin() { return bp->getYRangeMin(); }
double getYRangeMax() { return bp->getYRangeMax(); }
double getZRangeMin() { return bp->getZRangeMin(); }
double getZRangeMax() { return bp->getZRangeMax(); }
double getXRangeMin(unsigned short level) { return bp->getXRangeMin(); }
double getXRangeMax(unsigned short level) { return bp->getXRangeMax(); }
double getYRangeMin(unsigned short level) { return bp->getYRangeMin(); }
double getYRangeMax(unsigned short level) { return bp->getYRangeMax(); }
double getZRangeMin(unsigned short level) { return bp->getZRangeMin(); }
double getZRangeMax(unsigned short level) { return bp->getZRangeMax(); }
void test(PartBunch &bunch) { }
/// useful load balance information
void printLoadBalanceStats();
......
......@@ -60,12 +60,12 @@ public:
void applyConstantFocusing(PartBunch &bunch, double f, double r);
void test(PartBunch &bunch);
double getXRangeMin() {return 1.0;}
double getXRangeMax() {return 1.0;}
double getYRangeMin() {return 1.0;}
double getYRangeMax() {return 1.0;}
double getZRangeMin() {return 1.0;}
double getZRangeMax() {return 1.0;}
double getXRangeMin(unsigned short level) {return 1.0;}
double getXRangeMax(unsigned short level) {return 1.0;}
double getYRangeMin(unsigned short level) {return 1.0;}
double getYRangeMax(unsigned short level) {return 1.0;}
double getZRangeMin(unsigned short level) {return 1.0;}
double getZRangeMax(unsigned short level) {return 1.0;}
void computeAvgSpaceChargeForces(PartBunch &bunch);
void compute_temperature(PartBunch &bunch);
......
......@@ -4,6 +4,10 @@
//////////////////////////////////////////////////////////////
#include "Algorithms/PBunchDefs.h"
#ifdef HAVE_AMR_SOLVER
#include "Utilities/OpalException.h"
#endif
//////////////////////////////////////////////////////////////
class PartBunch;
//use Barton and Nackman Trick to avoid virtual functions
......@@ -21,14 +25,30 @@ public:
// given a charge-density field rho and a set of mesh spacings hr,
// compute the scalar potential in open space
virtual void computePotential(Field_t &rho, Vector_t hr) = 0;
/*!
* @param rho specifies the charge-density field
* @param baseLevel of adaptive mesh refinement solvers (AMR). Used in case of sub-cycling in time.
* @param finestLevel of AMR.
*/
#ifdef HAVE_AMR_SOLVER
virtual void solve(AmrFieldContainer_t &rho,
AmrFieldContainer_t &efield,
unsigned short baseLevel,
unsigned short finestLevel)
{
throw OpalException("PoissonSolver", "Not implemented.");
};
#endif
virtual void computePotential(Field_t &rho, Vector_t hr, double zshift) = 0;
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 double getXRangeMin(unsigned short level = 0) = 0;
virtual double getXRangeMax(unsigned short level = 0) = 0;
virtual double getYRangeMin(unsigned short level = 0) = 0;
virtual double getYRangeMax(unsigned short level = 0) = 0;
virtual double getZRangeMin(unsigned short level = 0) = 0;
virtual double getZRangeMax(unsigned short level = 0) = 0;
virtual void test(PartBunch &bunch) = 0 ;
virtual ~PoissonSolver(){};
......
Markdown is supported
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