3#include "dg/algorithm.h"
5#include "dg/backend/timer.h"
12template <
class ProductMPIGeometry,
class MIMatrix,
class LocalContainer>
13struct Fieldaligned< ProductMPIGeometry, MIMatrix, MPI_Vector<LocalContainer> >
16 template <
class Limiter>
18 const ProductMPIGeometry&
grid,
23 unsigned mx=12,
unsigned my=12,
24 double deltaPhi = -1, std::string interpolation_method =
"linear-nearest",
25 bool benchmark =
true):
27 mx, my,
deltaPhi, interpolation_method, benchmark)
30 template <
class Limiter>
32 const ProductMPIGeometry&
grid,
37 unsigned mx=12,
unsigned my=12,
38 double deltaPhi = -1, std::string interpolation_method =
"linear-nearest",
39 bool benchmark =
true);
40 template<
class ...Params>
61 void set_boundaries(
dg::bc bcz,
const MPI_Vector<LocalContainer>& left,
const MPI_Vector<LocalContainer>& right)
68 void set_boundaries(
dg::bc bcz,
const MPI_Vector<LocalContainer>& global,
double scal_left,
double scal_right)
76 void operator()(
enum whichMatrix which,
const MPI_Vector<LocalContainer>& in, MPI_Vector<LocalContainer>& out);
78 double deltaPhi()
const{
return m_deltaPhi;}
79 const MPI_Vector<LocalContainer>&
hbm()
const {
82 const MPI_Vector<LocalContainer>&
hbp()
const {
85 const MPI_Vector<LocalContainer>&
sqrtG()
const {
88 const MPI_Vector<LocalContainer>&
sqrtGm()
const {
91 const MPI_Vector<LocalContainer>&
sqrtGp()
const {
94 const MPI_Vector<LocalContainer>&
bphi()
const {
97 const MPI_Vector<LocalContainer>&
bphiM()
const {
100 const MPI_Vector<LocalContainer>&
bphiP()
const {
103 const MPI_Vector<LocalContainer>&
bbm()
const {
106 const MPI_Vector<LocalContainer>&
bbo()
const {
109 const MPI_Vector<LocalContainer>&
bbp()
const {
112 const ProductMPIGeometry&
grid()
const{
return *m_g;}
114 template<
class BinaryOp,
class UnaryOp>
115 MPI_Vector<LocalContainer>
evaluate( BinaryOp f, UnaryOp g,
unsigned p0,
116 unsigned rounds)
const;
117 std::string
method()
const{
return m_interpolation_method;}
119 void ePlus(
enum whichMatrix which,
const MPI_Vector<LocalContainer>& in, MPI_Vector<LocalContainer>& out);
120 void eMinus(
enum whichMatrix which,
const MPI_Vector<LocalContainer>& in, MPI_Vector<LocalContainer>& out);
121 void zero(
enum whichMatrix which,
const MPI_Vector<LocalContainer>& in, MPI_Vector<LocalContainer>& out);
122 MIMatrix m_plus, m_zero, m_minus, m_plusT, m_minusT;
123 MPI_Vector<LocalContainer> m_hbm, m_hbp;
124 MPI_Vector<LocalContainer> m_G, m_Gm, m_Gp;
125 MPI_Vector<LocalContainer> m_bphi, m_bphiM, m_bphiP;
126 MPI_Vector<LocalContainer> m_bbm, m_bbp, m_bbo;
128 MPI_Vector<LocalContainer> m_left, m_right;
129 MPI_Vector<LocalContainer> m_limiter;
130 mutable MPI_Vector<LocalContainer> m_ghostM, m_ghostP;
131 mutable std::vector<MPI_Vector<dg::View<const LocalContainer>> > m_f;
132 mutable std::vector<MPI_Vector<dg::View<LocalContainer>> > m_temp;
134 dg::bc m_bcx, m_bcy, m_bcz;
135 unsigned m_Nz, m_perp_size, m_mx, m_my;
136 double m_deltaPhi, m_eps;
137 std::string m_interpolation_method;
139 unsigned m_coords2, m_sizeZ;
141 mutable thrust::host_vector<double> m_buffer;
142 dg::detail::MPIContiguousGather m_from_minus, m_from_plus;
143 bool m_have_adjoint =
false;
144 void updateAdjoint( )
151 std::array<thrust::host_vector<double>,3> yp_trafo, ym_trafo;
152 thrust::host_vector<double>
hbp,
hbm;
153 thrust::host_vector<bool> in_boxp, in_boxm;
155 make_matrices( m_vec, grid_transform, global_grid_magnetic,
156 m_bcx, m_bcy, m_eps, m_mx, m_my, m_deltaPhi,
157 m_interpolation_method,
158 false,
true, vol2d0,
hbp,
hbm,
168 unsigned mx,
unsigned my,
169 double deltaPhi, std::string interpolation_method,
170 bool benchmark,
bool make_adjoint,
171 const MPI_Vector<thrust::host_vector<double>>& vol2d0,
172 thrust::host_vector<double>&
hbp,
173 thrust::host_vector<double>&
hbm,
174 thrust::host_vector<bool>& in_boxp,
175 thrust::host_vector<bool>& in_boxm,
176 std::array<thrust::host_vector<double>,3>& yp_trafo,
177 std::array<thrust::host_vector<double>,3>& ym_trafo
181 MPI_Comm_rank( m_g->communicator(), &rank);
182 std::string inter_m, project_m, fine_m;
183 detail::parse_method( interpolation_method, inter_m, project_m, fine_m);
184 if( benchmark && rank==0)
185 std::cout <<
"# Interpolation method: \""<<inter_m
186 <<
"\" projection method: \""<<project_m
187 <<
"\" fine grid \""<<fine_m<<
"\"\n";
191 if( benchmark) t.
tic();
197 grid_equidist_global.set( 1, grid_equidist_global.shape(0), grid_equidist_global.shape(1));
199 grid_magnetic->set( grid_transform->n() < 3 ? 4 : 7, grid_magnetic->Nx(), grid_magnetic->Ny());
200 global_grid_magnetic = grid_magnetic->global_geometry();
202 if( project_m !=
"dg" && fine_m ==
"dg")
204 unsigned rx = mx % grid_transform->nx(), ry = my % grid_transform->ny();
205 if( 0 != rx || 0 != ry)
207 if(rank==0)std::cerr <<
"#Warning: for projection method \"const\" mx and my must be multiples of nx and ny! Rounding up for you ...\n";
208 mx = mx + grid_transform->nx() - rx;
209 my = my + grid_transform->ny() - ry;
212 if( fine_m ==
"equi")
213 grid_fine_local.set( 1, grid_fine_local.shape(0), grid_fine_local.shape(1));
214 grid_fine_local.multiplyCellNumbers((
double)mx, (
double)my);
218 if(rank==0) std::cout <<
"# DS: High order grid gen took: "<<t.
diff()<<
"\n";
222 std::array<thrust::host_vector<double>,3> yp, ym;
223 detail::integrate_all_fieldlines2d( vec, *global_grid_magnetic,
224 grid_transform->local(), yp_trafo, vol2d0.data(),
hbp, in_boxp,
226 detail::integrate_all_fieldlines2d( vec, *global_grid_magnetic,
227 grid_transform->local(), ym_trafo, vol2d0.data(),
hbm, in_boxm,
233 grid_transform->local(),
dg::NEU,
dg::NEU, grid_transform->n() < 3 ?
"cubic" :
"dg");
236 for(
int i=0; i<2; i++)
245 if(rank==0) std::cout <<
"# DS: Fieldline integration took: "<<t.
diff()<<
"\n";
251 if( project_m ==
"dg")
256 std::array<dg::HVec*,3> xcomp{ &yp[0], &Xf, &ym[0]};
257 std::array<dg::HVec*,3> ycomp{ &yp[1], &Yf, &ym[1]};
258 std::array<MIMatrix*,3> result{ &m_plus, &m_zero, &m_minus};
259 std::array<MIMatrix*,3> resultT{ &m_plusT, &m_zero, &m_minusT};
260 for(
unsigned u=0; u<3; u++)
265 grid_transform->global(),
bcx,
bcy,
"dg");
271 grid_equidist_global,
bcx,
bcy, inter_m) *
276 if( project_m !=
"dg")
282 if( make_adjoint and u != 1)
293 m_have_adjoint =
true;
300 if(rank==0) std::cout <<
"# DS: Multiplication PI took: "<<t.
diff()<<
"\n";
305template<
class MPIGeometry,
class MIMatrix,
class LocalContainer>
306template <
class Limiter>
307Fieldaligned<MPIGeometry, MIMatrix, MPI_Vector<LocalContainer> >::Fieldaligned(
309 const MPIGeometry& grid,
311 unsigned mx,
unsigned my,
312 double deltaPhi, std::string interpolation_method,
bool benchmark
314 m_g(grid), m_bcx(bcx), m_bcy(bcy), m_bcz(grid.bcz()),
315 m_Nz( grid.local().Nz()), m_mx(mx), m_my(my), m_eps(eps),
316 m_interpolation_method(interpolation_method),
320 MPI_Comm_rank( grid.communicator(), &rank);
321 int dims[3], periods[3], coords[3];
322 MPI_Cart_get( m_g->communicator(), 3, dims, periods, coords);
323 m_coords2 = coords[2], m_sizeZ = dims[2];
326 if( (grid.bcx() == PER && bcx != PER) || (grid.bcx() != PER && bcx == PER) )
327 throw(
dg::Error(
dg::Message(_ping_)<<
"Fieldaligned: Got conflicting periodicity in x. The grid says "<<
bc2str(grid.bcx())<<
" while the parameter says "<<
bc2str(bcx)));
328 if( (grid.bcy() == PER && bcy != PER) || (grid.bcy() != PER && bcy == PER) )
329 throw(
dg::Error(
dg::Message(_ping_)<<
"Fieldaligned: Got conflicting boundary conditions in y. The grid says "<<
bc2str(grid.bcy())<<
" while the parameter says "<<
bc2str(bcy)));
330 if( deltaPhi <=0) deltaPhi = grid.hz();
331 m_deltaPhi = deltaPhi;
338 std::array<thrust::host_vector<double>,3> yp_trafo, ym_trafo;
339 thrust::host_vector<double> hbp, hbm;
340 thrust::host_vector<bool> in_boxp, in_boxm;
342 make_matrices( vec, grid_transform, global_grid_magnetic,
343 bcx, bcy, eps, mx, my, m_deltaPhi, interpolation_method,
344 benchmark,
false, vol2d0, hbp, hbm,
348 dg::HVec hbphi( yp_trafo[2]), hbphiP(hbphi), hbphiM(hbphi);
354 for(
unsigned i=0; i<hbphiP.size(); i++)
356 hbphiP[i] = vec.
z()(yp_trafo[0][i], yp_trafo[1][i]);
357 hbphiM[i] = vec.
z()(ym_trafo[0][i], ym_trafo[1][i]);
364 for(
unsigned i=0; i<yp_trafo[0].size(); i++)
367 yp_trafo[1][i], *global_grid_magnetic);
369 ym_trafo[1][i], *global_grid_magnetic);
386 m_f =
dg::split( (
const MPI_Vector<LocalContainer>&)m_hbm, grid);
391 thrust::host_vector<double> bbm( in_boxp.size(),0.), bbo(bbm), bbp(bbm);
392 for(
unsigned i=0; i<in_boxp.size(); i++)
394 if( !in_boxp[i] && !in_boxm[i])
396 else if( !in_boxp[i] && in_boxm[i])
398 else if( in_boxp[i] && !in_boxm[i])
408 m_perp_size = grid_transform->local().size();
411 m_ghostM = m_ghostP = m_right = m_left;
414 dg::detail::MsgChunk chunk { 0, (int)m_perp_size};
416 MPI_Cart_shift( grid.comm(2), 0, +1, &source, &dest);
417 std::map<int, thrust::host_vector<dg::detail::MsgChunk>> recvMsgP =
418 {{ dest, thrust::host_vector<dg::detail::MsgChunk>( 1, chunk)}};
419 m_from_plus = dg::detail::MPIContiguousGather( recvMsgP, grid.comm(2));
421 MPI_Cart_shift( grid.comm(2), 0, -1, &source, &dest);
422 std::map<int, thrust::host_vector<dg::detail::MsgChunk>> recvMsgM =
423 {{ dest, thrust::host_vector<dg::detail::MsgChunk>( 1, chunk)}};
424 m_from_minus = dg::detail::MPIContiguousGather( recvMsgM, grid.comm(2));
428template<
class G,
class M,
class container>
429void Fieldaligned<G, M, MPI_Vector<container> >::operator()(
enum
431 MPI_Vector<container>& fe)
433 if(which == einsPlus || which == einsMinusT) ePlus( which, f, fe);
434 if(which == einsMinus || which == einsPlusT) eMinus( which, f, fe);
435 if( which == einsPlus || which == einsMinusT ) ePlus( which, f, fe);
436 else if(which == einsMinus || which == einsPlusT ) eMinus( which, f, fe);
437 else if(which == zeroMinus || which == zeroPlus ||
438 which == zeroMinusT|| which == zeroPlusT ||
439 which == zeroForw )
zero( which, f, fe);
441template<
class G,
class M,
class container>
442void Fieldaligned<G, M, MPI_Vector<container> >::zero(
enum whichMatrix which,
const MPI_Vector<container>& f, MPI_Vector<container>& f0)
447 for(
unsigned i0=0; i0<m_Nz; i0++)
449 if(which == zeroPlus)
451 else if(which == zeroMinus)
453 else if(which == zeroPlusT)
455 if( ! m_have_adjoint) updateAdjoint( );
458 else if(which == zeroMinusT)
460 if( ! m_have_adjoint) updateAdjoint( );
463 else if( which == zeroForw)
465 if ( m_interpolation_method !=
"dg" )
475template<
class G,
class M,
class container>
476void Fieldaligned<G,M, MPI_Vector<container> >::ePlus(
enum whichMatrix which,
const MPI_Vector<container>& f, MPI_Vector<container>& fpe )
480 MPI_Vector<dg::View<container>> send_buf(
481 {m_ghostP.data().data(), m_ghostP.size()}, m_g->get_perp_comm());
483 for(
unsigned i0=0; i0<m_Nz; i0++)
486 bool comm_plane = (m_sizeZ != 1 and i0 == m_Nz -1);
487 unsigned ip = (i0==m_Nz-1) ? 0:i0+1;
488 if(which == einsPlus)
490 else if(which == einsMinusT)
492 if( ! m_have_adjoint) updateAdjoint( );
493 dg::blas2::symv( m_minusT, m_f[ip], comm_plane ? send_buf : m_temp[i0]);
500 unsigned i0 = m_Nz-1;
501 m_from_plus.global_gather_init( send_buf.data(), m_temp[i0].data());
502 m_from_plus.global_gather_wait( m_temp[i0].data());
507 if( m_bcz !=
dg::PER && m_g->local().z1() == m_g->global().z1())
519template<
class G,
class M,
class container>
520void Fieldaligned<G, M, MPI_Vector<container> >::eMinus(
enum
521 whichMatrix which,
const MPI_Vector<container>& f, MPI_Vector<container>& fme )
524 MPI_Comm_rank(m_g->communicator(), &rank);
527 MPI_Vector<dg::View<container>> send_buf(
528 {m_ghostM.data().data(), m_ghostM.size()}, m_g->get_perp_comm());
530 for(
unsigned i0=0; i0<m_Nz; i0++)
533 bool comm_plane = (m_sizeZ != 1 and i0 == 0);
534 unsigned im = (i0==0) ? m_Nz-1:i0-1;
535 if(which == einsPlusT)
537 if( ! m_have_adjoint) updateAdjoint( );
538 dg::blas2::symv( m_plusT, m_f[im], comm_plane ? send_buf : m_temp[i0]);
540 else if(which == einsMinus)
541 dg::blas2::symv( m_minus, m_f[im], comm_plane ? send_buf : m_temp[i0]);
548 m_from_minus.global_gather_init( send_buf.data(), m_temp[i0].data());
549 m_from_minus.global_gather_wait( m_temp[i0].data());
554 if( m_bcz !=
dg::PER && m_g->local().z0() == m_g->global().z0())
566template<
class G,
class M,
class container>
567template<
class BinaryOp,
class UnaryOp>
568MPI_Vector<container> Fieldaligned<G,M, MPI_Vector<container> >::evaluate( BinaryOp binary, UnaryOp unary,
unsigned p0,
unsigned rounds)
const
572 assert( p0 < m_g->global().Nz());
574 MPI_Vector<container> init2d =
dg::pullback( binary, *g2d);
576 unsigned globalNz = m_g->global().Nz();
578 MPI_Vector<container> temp(init2d), tempP(init2d), tempM(init2d);
580 std::vector<MPI_Vector<container> > plus2d(globalNz, zero2d), minus2d(plus2d), result(plus2d);
581 unsigned turns = rounds;
582 if( turns ==0) turns++;
584 for(
unsigned r=0; r<turns; r++)
585 for(
unsigned i0=0; i0<globalNz; i0++)
589 unsigned rep = r*globalNz + i0;
590 for(
unsigned k=0; k<rep; k++)
607 for(
unsigned i0=0; i0<m_Nz; i0++)
609 int idx = (int)(i0+m_coords2*m_Nz) - (int)p0;
611 result[i0] = plus2d[idx];
613 result[i0] = minus2d[abs(idx)];
614 thrust::copy( result[i0].data().begin(), result[i0].data().end(),
615 vec3d.data().begin() + i0*m_perp_size);
620 for(
unsigned i0=0; i0<globalNz; i0++)
622 unsigned revi0 = (globalNz - i0)%globalNz;
627 for(
unsigned i0=0; i0<m_Nz; i0++)
629 int idx = ((int)i0 + m_coords2*m_Nz -(int)p0 + globalNz)%globalNz;
630 thrust::copy( result[idx].data().begin(), result[idx].data().end(),
631 vec3d.data().begin() + i0*m_perp_size);
644template<
class BinaryOp,
class UnaryOp>
648 const BinaryOp& binary,
649 const UnaryOp& unary,
654 unsigned Nz = grid.
Nz();
659 std::vector<dg::MHVec> plus2d(Nz, tempP), minus2d(plus2d), result(plus2d);
661 std::array<dg::HVec,3> yy0{
665 dg::geo::detail::DSFieldCylindrical3 cyl_field(vec);
666 double deltaPhi = grid.
hz();
667 double phiM0 = 0., phiP0 = 0.;
668 unsigned turns = rounds;
669 if( turns == 0) turns++;
670 for(
unsigned r=0; r<turns; r++)
671 for(
unsigned i0=0; i0<Nz; i0++)
673 unsigned rep = r*Nz + i0;
675 tempM = tempP = init2d;
679 "Dormand-Prince-7-4-5", std::array<double,3>{0,0,0});
682 for(
unsigned i=0; i<g2d->local().size(); i++)
685 double phiM1 = phiM0 + deltaPhi;
687 coords0{yy0[0][i],yy0[1][i],yy0[2][i]}, coords1;
688 odeint.integrate_in_domain( phiM0, coords0, phiM1,
689 coords1, deltaPhi, g2d->global(), eps);
690 yy1[0][i] = coords1[0], yy1[1][i] = coords1[1], yy1[2][i] =
692 tempM.data()[i] = binary( yy1[0][i], yy1[1][i]);
695 double phiP1 = phiP0 - deltaPhi;
696 coords0 = std::array<double,3>{xx0[0][i],xx0[1][i],xx0[2][i]};
697 odeint.integrate_in_domain( phiP0, coords0, phiP1,
698 coords1, -deltaPhi, g2d->global(), eps);
699 xx1[0][i] = coords1[0], xx1[1][i] = coords1[1], xx1[2][i] =
701 tempP.
data()[i] = binary( xx1[0][i], xx1[1][i]);
703 std::swap( yy0, yy1);
704 std::swap( xx0, xx1);
714 int dims[3], periods[3], coords[3];
715 MPI_Cart_get( grid.
communicator(), 3, dims, periods, coords);
716 unsigned coords2 = coords[2];
719 for(
unsigned i0=0; i0<grid.
local().Nz(); i0++)
721 int idx = (int)(i0+coords2*grid.
local().Nz()) - (int)p0;
723 result[i0] = plus2d[idx];
725 result[i0] = minus2d[abs(idx)];
726 thrust::copy( result[i0].data().begin(), result[i0].data().end(),
727 vec3d.
data().begin() + i0*g2d->local().size());
732 for(
unsigned i0=0; i0<Nz; i0++)
734 unsigned revi0 = (Nz - i0)%Nz;
739 for(
unsigned i0=0; i0<grid.
local().Nz(); i0++)
741 int idx = ((int)i0 +coords2*grid.
local().Nz()-(
int)p0 + Nz)%Nz;
743 thrust::copy( result[idx].data().begin(), result[idx].data().end(),
744 vec3d.
data().begin() + i0*g2d->local().size());
DG_DEVICE T zero(T x, Ts ...xs)
DG_DEVICE double cooY2d(double x, double y)
DG_DEVICE double cooX2d(double x, double y)
void copy(const ContainerTypeIn &source, ContainerTypeOut &target)
void axpby(value_type alpha, const ContainerType1 &x, value_type1 beta, ContainerType &y)
void pointwiseDot(value_type alpha, const ContainerType1 &x1, const ContainerType2 &x2, value_type1 beta, ContainerType &y)
void assign(const from_ContainerType &from, ContainerType &to, Params &&... ps)
void scal(ContainerType &x, value_type alpha)
void transfer(const MatrixType &x, AnotherMatrixType &y)
void symv(MatrixType &&M, const ContainerType1 &x, ContainerType2 &y)
std::string bc2str(bc bcx)
auto weights(const Topology &g)
auto evaluate(Functor &&f, const Topology &g)
whichMatrix
Enum for the use in Fieldaligned.
Definition fieldaligned.h:16
thrust::host_vector< double > fieldaligned_evaluate(const aProductGeometry3d &grid, const CylindricalVectorLvl0 &vec, const BinaryOp &binary, const UnaryOp &unary, unsigned p0, unsigned rounds, double eps=1e-5)
Evaluate a 2d functor and transform to all planes along the fieldlines
Definition fieldaligned.h:1014
ONE FullLimiter
Full Limiter means there is a limiter everywhere.
Definition fieldaligned.h:30
Topology::host_vector forward_transform(const typename Topology::host_vector &in, const Topology &g)
real_type interpolate(dg::space sp, const host_vector &v, real_type x, const RealGrid1d< real_type > &g, dg::bc bcx=dg::NEU)
dg::MIHMatrix_t< typename MPITopology::value_type > projection(const MPITopology &g_new, const MPITopology &g_old, std::string method="dg")
dg::SparseMatrix< int, real_type, thrust::host_vector > interpolation(const RecursiveHostVector &x, const aRealTopology< real_type, Nd > &g, std::array< dg::bc, Nd > bcx, std::string method="dg")
CylindricalVectorLvl1 createBHat(const TokamakMagneticField &mag)
Contravariant components of the magnetic unit vector field and its Divergence and derivative in cylin...
Definition magnetic_field.h:934
void convertLocal2GlobalCols(dg::IHMatrix_t< real_type > &local, const ConversionPolicy &policy)
dg::IHMatrix_t< real_type > convertGlobal2LocalRows(const dg::IHMatrix_t< real_type > &global, const ConversionPolicy &row_policy)
dg::MIHMatrix_t< real_type > make_mpi_matrix(const dg::IHMatrix_t< real_type > &global_cols, const ConversionPolicy &col_policy)
Geometry::host_vector pullback(const Functor &f, const Geometry &g)
void assign3dfrom2d(const host_vector &in2d, Container &out, const Topology &grid)
dg::IHMatrix_t< real_type > inv_backproject(const aRealTopology< real_type, Nd > &g)
dg::IHMatrix_t< real_type > backproject(const aRealTopology< real_type, Nd > &g)
void split(SharedContainer &in, std::vector< View< SharedContainer > > &out, const aRealTopology3d< real_type > &grid)
ContainerType volume(const SparseTensor< ContainerType > &t)
thrust::host_vector< double > HVec
const container & data() const
SparseMatrix transpose() const
MPI_Comm communicator() const
const RealGrid< real_type, Nd > & local() const
aRealMPIGeometry2d< real_type > * perp_grid() const
Definition fluxfunctions.h:416
This struct bundles a vector field and its divergence.
Definition fluxfunctions.h:444
const CylindricalFunctor & z() const
z-component of the vector
Definition fluxfunctions.h:474
dg::bc bcx() const
Definition fieldaligned.h:322
const container & hbp() const
Distance between the planes .
Definition fieldaligned.h:399
const container & bphi() const
bphi
Definition fieldaligned.h:419
const container & bbp() const
Mask plus, 1 if fieldline intersects wall in plus direction but not in minus direction,...
Definition fieldaligned.h:444
container evaluate(BinaryOp binary, UnaryOp unary, unsigned p0, unsigned rounds) const
Evaluate a 2d functor and transform to all planes along the fieldline
const container & sqrtG() const
Volume form (including weights) .
Definition fieldaligned.h:404
void operator()(enum whichMatrix which, const container &in, container &out)
Apply the interpolation to three-dimensional vectors.
const container & sqrtGm() const
Volume form on minus plane (including weights) .
Definition fieldaligned.h:409
void set_boundaries(dg::bc bcz, double left, double right)
Set boundary conditions in the limiter region.
Definition fieldaligned.h:339
std::string method() const
Return the interpolation_method string given in the constructor.
Definition fieldaligned.h:509
const container & hbm() const
Distance between the planes and the boundary .
Definition fieldaligned.h:394
const container & bphiM() const
bphi on minus plane
Definition fieldaligned.h:424
dg::bc bcy() const
Definition fieldaligned.h:325
double deltaPhi() const
Definition fieldaligned.h:391
Fieldaligned()
do not allocate memory; no member call except construct is valid
Definition fieldaligned.h:274
void construct(Params &&...ps)
Perfect forward parameters to one of the constructors.
Definition fieldaligned.h:316
const container & bbo() const
Mask both, 1 if fieldline intersects wall in plus direction and in minus direction,...
Definition fieldaligned.h:439
const container & bbm() const
Mask minus, 1 if fieldline intersects wall in minus direction but not in plus direction,...
Definition fieldaligned.h:434
const container & sqrtGp() const
Volume form on plus plane (including weights) .
Definition fieldaligned.h:414
const container & bphiP() const
bphi on plus plane
Definition fieldaligned.h:429
const ProductGeometry & grid() const
Grid used for construction.
Definition fieldaligned.h:448
A tokamak field as given by R0, Psi and Ipol plus Meta-data like shape and equilibrium.
Definition magnetic_field.h:165