3#include "dg/algorithm.h"
69 DSCentered(
double alpha,
double beta) : m_alpha(alpha), m_beta(beta){}
71 void operator()(
double& dsf,
double fm,
double fo,
double fp,
double hm,
75 fm*( 1./(hp+hm) - 1./hm) +
77 fp*( 1./hp - 1./(hp+hm))
82 double m_alpha, m_beta;
86 DSSCentered(
double alpha,
double beta,
double delta) : m_alpha(alpha), m_beta(beta), m_delta(delta){}
88 void operator()(
double& dssf,
double fm,
double fo,
double fp,
double hm,
92 2.*fm/(hp+hm)/hm - 2.*fo/hp/hm + 2.*fp/(hp+hm)/hp
97 void operator()(
double& dssf,
double fm,
double fo,
double fp,
98 double bPm,
double bP0,
double bPp)
100 double bP2 = (bPp+bP0)/2.;
101 double bM2 = (bPm+bP0)/2.;
102 double fm2 = (fo-fm)/m_delta;
103 double fp2 = (fp-fo)/m_delta;
105 dssf = m_alpha*bP0*( bP2*fp2 - bM2*fm2)/m_delta + m_beta*dssf;
109 double m_alpha, m_beta, m_delta;
113 DSSDCentered(
double alpha,
double beta,
double delta) :
114 m_alpha( alpha), m_beta(beta), m_delta(delta){}
116 void operator()(
double& dssdf,
double fm,
double fo,
double fp,
double Gm,
117 double Go,
double Gp,
double bPm,
double bP0,
double bPp)
120 double bP2 = (bPp+bP0)/2.;
121 double bM2 = (bPm+bP0)/2.;
122 double fm2 = (fo-fm)/m_delta;
123 double fp2 = (fp-fo)/m_delta;
124 double gp2 = (Gp + Go)/Go/2.;
125 double gm2 = (Gm + Go)/Go/2.;
127 dssdf = m_alpha*( gp2*fp2*bP2*bP2 - bM2*bM2*gm2*fm2)/m_delta + m_beta*dssdf;
135 double m_alpha, m_beta, m_delta;
169template<
class FieldAligned,
class container>
171 const container& f,
const container& fp, container& fmg, container& fpg,
172 dg::bc bound, std::array<double,2> boundary_value = {0,0})
174 double delta = fa.deltaPhi();
177 double dbm = boundary_value[0], dbp = boundary_value[1];
179 double fp,
double& fmg,
double& fpg,
180 double hbm,
double hbp,
double bbm,
double bbo,
double bbp
183 double hm = delta, hp = delta;
184 double plus=0, minus=0, bothP=0, bothM = 0;
185 plus = dbp*hp*(hm+hp)/(2.*hbp+hm) +
186 fo*(2.*hbp+hm-hp)*(hm+hp)/hm/(2.*hbp+hm) + fm*hp*(-2.*hbp +
187 hp)/hm/(2.*hbp + hm);
188 minus = fp*hm*(-2.*hbm+hm)/hp/(2.*hbm+hp) -
189 dbm*hm*(hm+hp)/(2.*hbm+hp) +
190 fo*(2.*hbm-hm+hp)*(hm+hp)/hp/(2.*hbm+hp);
191 bothM = fo + dbp*hm*(-2.*hbm + hm)/2./(hbm+hbp) -
192 dbm*hm*(2.*hbp+hm)/2./(hbm+hbp);
193 bothP = fo + dbp*hp*(2.*hbm + hp)/2./(hbm+hbp) +
194 dbm*hp*(2.*hbp-hp)/2./(hbm+hbp);
195 fmg = (1.-bbo-bbm)*fm + bbm*minus + bbo*bothM;
196 fpg = (1.-bbo-bbp)*fp + bbp*plus + bbo*bothP;
197 }, fm, f, fp, fmg, fpg, fa.hbm(), fa.hbp(), fa.bbm(),
198 fa.bbo(), fa.bbp() );
202 double fbm = boundary_value[0], fbp = boundary_value[1];
204 double fp,
double& fmg,
double& fpg,
205 double hbm,
double hbp,
double bbm,
double bbo,
double bbp
208 double hm = delta, hp = delta;
209 double plus=0, minus=0, bothP=0, bothM = 0;
210 plus = fm*hp*(-hbp + hp)/hm/(hbp+hm) + fo*(hbp-hp)*(hm+hp)/hbp/hm
211 +fbp*hp*(hm+hp)/hbp/(hbp+hm);
212 minus = +fo*(hbm-hm)*(hm+hp)/hbm/hp + fbm*hm*(hm+hp)/hbm/(hbm+hp)
213 + fp*hm*(-hbm+hm)/hp/(hbm+hp);
214 bothM = fbp*hm*(-hbm+hm)/hbp/(hbm+hbp) +
215 fo*(hbm-hm)*(hbp+hm)/hbm/hbp + fbm*hm*(hbp+hm)/hbm/(hbm+hbp);
216 bothP = fo*(hbp-hp)*(hbm+hp)/hbm/hbp +
217 fbp*hp*(hbm+hp)/hbp/(hbm+hbp) + fbm*hp*(-hbp+hp)/hbm/(hbm+hbp);
218 fmg = (1.-bbo-bbm)*fm + bbm*minus + bbo*bothM;
219 fpg = (1.-bbo-bbp)*fp + bbp*plus + bbo*bothP;
220 }, fm, f, fp, fmg, fpg, fa.hbm(), fa.hbp(), fa.bbm(),
241template<
class FieldAligned,
class container>
243 const container& fp, container& fmg, container& fpg,
244 dg::bc bound, std::array<double,2> boundary_value = {0,0})
246 double delta = fa.deltaPhi();
249 double dbm = boundary_value[0], dbp = boundary_value[1];
251 double& fmg,
double& fpg,
double bbm,
double bbp
253 double hm = delta, hp = delta;
254 double plus=0, minus=0;
255 plus = fm + dbp*(hp+hm);
256 minus = fp - dbm*(hp+hm);
257 fmg = (1.-bbm)*fm + bbm*minus;
258 fpg = (1.-bbp)*fp + bbp*plus;
259 }, fm, fp, fmg, fpg, fa.bbm(), fa.bbp() );
263 double fbm = boundary_value[0], fbp = boundary_value[1];
265 double& fmg,
double& fpg,
double hbm,
266 double hbp,
double bbm,
double bbo,
double bbp
268 double hm = delta, hp = delta;
269 double plus=0, minus=0, bothP=0, bothM = 0;
270 plus = fm + (fbp-fm)/(hbp+hm)*(hp+hm) ;
271 minus = fp - (hp+hm)*(fp-fbm)/(hp+hbm);
272 bothM = fbp + (fbp-fbm)/(hbp+hbm)*(hp+hbm);
273 bothP = fbp - (fbp-fbm)/(hbp+hbm)*(hbp+hm);
274 fmg = (1.-bbo-bbm)*fm + bbm*minus + bbo*bothM;
275 fpg = (1.-bbo-bbp)*fp + bbp*plus + bbo*bothP;
276 }, fm, fp, fmg, fpg, fa.hbm(), fa.hbp(), fa.bbm(),
306template<
class FieldAligned,
class container>
308 const container& fp, container& fmg, container& fpg)
311 double& fmg,
double& fpg,
312 double bbm,
double bbo,
double bbp
314 fmg = (1.-bbo-bbm)*fm + (bbm+bbo)*(-fm);
315 fpg = (1.-bbo-bbp)*fp + (bbp+bbo)*(-fp);
316 }, fm, fp, fmg, fpg, fa.bbm(), fa.bbo(), fa.bbp() );
347template<
class ProductGeometry,
class IMatrix,
class container >
361 template <
class Limiter>
367 unsigned mx=10,
unsigned my=10,
368 double deltaPhi=-1, std::string interpolation_method =
"dg",
369 bool benchmark=
true):
370 DS(
FA( vec, grid, bcx, bcy, limit, eps, mx, my, deltaPhi,
371 interpolation_method, benchmark) )
381 template<
class Limiter>
387 unsigned mx=10,
unsigned my=10,
388 double deltaPhi=-1, std::string interpolation_method =
"dg",
389 bool benchmark=
true):
390 DS(
FA( vec, grid, bcx, bcy, limit, eps, mx, my, deltaPhi,
391 interpolation_method, benchmark))
405 template<
class ...Params>
409 *
this =
DS( std::forward<Params>( ps)...);
414 m_fa.set_boundaries( bcz, left, right);
418 m_fa.set_boundaries( bcz, left, right);
422 m_fa.set_boundaries( bcz, global, scal_left, scal_right);
431 void forward(
double alpha,
const container& f,
double beta, container& g){
434 ds_forward( m_fa, alpha, m_tempO, m_tempP, beta, g);
442 void forward2(
double alpha,
const container& f,
double beta, container& g){
446 ds_forward2( m_fa, alpha, m_tempO, m_tempP, m_tempM, beta, g);
454 void backward(
double alpha,
const container& f,
double beta, container& g){
457 ds_backward( m_fa, alpha, m_tempM, m_tempO, beta, g);
465 void backward2(
double alpha,
const container& f,
double beta, container& g){
469 ds_backward2( m_fa, alpha, m_tempP, m_tempM, m_tempO, beta, g);
481 void centered(
double alpha,
const container& f,
double beta, container& g){
484 ds_centered( m_fa, alpha, m_tempM, m_tempP, beta, g);
488 double alpha,
const container& f,
double beta, container& g,
dg::bc bound,
489 std::array<double,2> boundary_value = {0,0}){
493 bound, boundary_value);
494 ds_centered( m_fa, alpha, m_tempM, m_tempP, beta, g);
512 void forward(
const container& f, container& g){
528 void divForward(
double alpha,
const container& f,
double beta, container& g){
536 void divBackward(
double alpha,
const container& f,
double beta, container& g){
544 void divCentered(
double alpha,
const container& f,
double beta, container& g){
575 ds(dir, 1., f, 0., g);
583 void ds(
dg::direction dir,
double alpha,
const container& f,
double beta, container& g);
592 div(dir, 1., f, 0., g);
612 void symv(
const container& f, container& g){
symv( 1., f, 0., g);}
622 void symv(
double alpha,
const container& f,
double beta, container& g);
629 void dss(
const container& f, container& g){
637 void dss(
double alpha,
const container& f,
double beta, container& g){
641 dss_centered( m_fa, alpha, m_tempM, m_tempO, m_tempP, beta, g);
645 double alpha,
const container& f,
double beta, container& g,
dg::bc bound,
646 std::array<double,2> boundary_value = {0,0}){
651 bound, boundary_value);
652 dss_centered( m_fa, alpha, m_tempM, m_tempO, m_tempP, beta, g);
655 void dssd(
double alpha,
const container& f,
double
660 dssd_centered( m_fa, alpha, m_tempM, m_tempO, m_tempP, beta, g);
664 container& f,
double beta, container& g,
dg::bc bound,
665 std::array<double,2> boundary_value = {0,0}){
670 bound, boundary_value);
688 container m_tempP, m_tempO, m_tempM;
694template<
class Geometry,
class I,
class container>
697 m_tempP = fa.sqrtG(), m_tempM = m_tempO = m_tempP;
700template<
class G,
class I,
class container>
702 const container& f,
double beta, container& dsf) {
705 return centered( alpha, f, beta, dsf);
707 return forward( alpha, f, beta, dsf);
709 return backward( alpha, f, beta, dsf);
712template<
class G,
class I,
class container>
714 const container& f,
double beta, container& dsf) {
717 return divCentered( alpha, f, beta, dsf);
719 return divForward( alpha, f, beta, dsf);
721 return divBackward( alpha, f, beta, dsf);
726template<
class G,
class I,
class container>
729 dssd( alpha, f, beta, dsTdsf);
743template<
class FieldAligned,
class container>
744void ds_forward(
const FieldAligned& fa,
double alpha,
const container& f,
745 const container& fp,
double beta, container& g)
748 double delta = fa.deltaPhi();
750 double& dsf,
double fo,
double fp,
double bphi){
751 dsf = alpha*bphi*( fp - fo)/delta + beta*dsf;
753 g, f, fp, fa.bphi());
765template<
class FieldAligned,
class container>
766void ds_forward2(
const FieldAligned& fa,
double alpha,
const container& f,
767 const container& fp,
const container& fpp,
double beta, container& g)
770 double delta = fa.deltaPhi();
772 double& dsf,
double fo,
double fp,
double fpp,
double bphi){
773 dsf = alpha*bphi*( -3.*fo + 4.*fp - fpp)/2./delta
776 g, f, fp, fpp, fa.bphi());
789template<
class FieldAligned,
class container>
790void ds_backward(
const FieldAligned& fa,
double alpha,
const container& fm,
791 const container& f,
double beta, container& g)
794 double delta = fa.deltaPhi();
796 double& dsf,
double fo,
double fm,
double bphi){
797 dsf = alpha*bphi*( fo - fm)/delta + beta*dsf;
799 g, f, fm, fa.bphi());
812template<
class FieldAligned,
class container>
813void ds_backward2(
const FieldAligned& fa,
double alpha,
const container& fmm,
const container& fm,
const container& f,
double beta, container& g)
816 double delta = fa.deltaPhi();
818 double& dsf,
double fo,
double fm,
double fmm,
double bphi){
819 dsf = alpha*bphi*( 3.*fo - 4.*fm + fmm)/2./delta
822 g, f, fm, fmm, fa.bphi());
841template<
class FieldAligned,
class container>
842void ds_centered(
const FieldAligned& fa,
double alpha,
const container& fm,
843 const container& fp,
double beta, container& g)
846 double delta=fa.deltaPhi();
848 double fp,
double bphi){
849 g = alpha*bphi*(fp-fm)/2./delta + beta*g;
850 }, g, fm, fp, fa.bphi());
863template<
class FieldAligned,
class container>
864void dss_centered(
const FieldAligned& fa,
double alpha,
const container& fm,
865 const container& f,
const container& fp,
double beta, container& g)
868 g, fm, f, fp, fa.bphiM(), fa.bphi(), fa.bphiP());
881template<
class FieldAligned,
class container>
882void dssd_centered(
const FieldAligned& fa,
double alpha,
const container& fm,
883 const container& f,
const container& fp,
double beta, container& g)
886 g, fm, f, fp, fa.sqrtGm(), fa.sqrtG(), fa.sqrtGp(),
887 fa.bphiM(), fa.bphi(), fa.bphiP());
903template<
class FieldAligned,
class container>
905 const container& f,
double beta, container& g)
907 double delta = fa.deltaPhi();
909 double f1,
double Gm,
double G0,
double bPm,
double bP0){
910 dsf = alpha*(bP0*G0*f0 - bPm*Gm*f1)/G0/delta + beta*dsf; },
911 g, f, fm, fa.sqrtGm(), fa.sqrtG(), fa.bphiM(), fa.bphi());
927template<
class FieldAligned,
class container>
928void ds_divForward(
const FieldAligned& fa,
double alpha,
const container& f,
929 const container& fp,
double beta, container& g)
931 double delta = fa.deltaPhi();
933 double f1,
double Gp,
double G0,
double bPp,
double bP0){
934 dsf = alpha*(bPp*Gp*f1 - bP0*G0*f0)/G0/delta + beta*dsf; },
935 g, f, fp, fa.sqrtGp(), fa.sqrtG(), fa.bphiP(), fa.bphi());
950template<
class FieldAligned,
class container>
951void ds_divCentered(
const FieldAligned& fa,
double alpha,
const container& fm,
const container& fp,
952 double beta, container& g)
954 double delta = fa.deltaPhi();
956 double fp,
double Gm,
double Gp,
double G0,
957 double bPm,
double bP0,
double bPp)
959 dsf = alpha*( fp*Gp*bPp - fm*Gm*bPm )/G0/2./delta + beta*dsf;
960 }, g, fm, fp, fa.sqrtGm(),
961 fa.sqrtGp(), fa.sqrtG(), fa.bphiM(), fa.bphi(), fa.bphiP());
980template<
class FieldAligned,
class container>
982 const container& fm,
const container& fp,
double beta, container& g)
986 g = alpha*(fp+fm)/2. + beta*g;
1005template<
class FieldAligned,
class container>
1007 const container& fm,
const container& fp,
double beta, container& g)
1016template<
class G,
class I,
class V>
1017struct TensorTraits< geo::DS<G,I, V> >
void plus(ContainerType &x, value_type alpha)
void subroutine(Subroutine f, ContainerType &&x, ContainerTypes &&... xs)
void dssd_centered(const FieldAligned &fa, double alpha, const container &fm, const container &f, const container &fp, double beta, container &g)
Centered derivative .
Definition ds.h:882
void ds_divForward(const FieldAligned &fa, double alpha, const container &f, const container &fp, double beta, container &g)
forward derivative
Definition ds.h:928
void assign_bc_along_field_1st(const FieldAligned &fa, const container &fm, const container &fp, container &fmg, container &fpg, dg::bc bound, std::array< double, 2 > boundary_value={0, 0})
Assign boundary conditions along magnetic field lines interpolating a 1st order polynomial (a line)
Definition ds.h:242
void assign_bc_along_field_2nd(const FieldAligned &fa, const container &fm, const container &f, const container &fp, container &fmg, container &fpg, dg::bc bound, std::array< double, 2 > boundary_value={0, 0})
Assign boundary conditions along magnetic field lines interpolating a 2nd order polynomial.
Definition ds.h:170
void ds_backward(const FieldAligned &fa, double alpha, const container &fm, const container &f, double beta, container &g)
backward derivative
Definition ds.h:790
void ds_centered(const FieldAligned &fa, double alpha, const container &fm, const container &fp, double beta, container &g)
centered derivative
Definition ds.h:842
void ds_divBackward(const FieldAligned &fa, double alpha, const container &fm, const container &f, double beta, container &g)
backward derivative
Definition ds.h:904
void ds_average(const FieldAligned &fa, double alpha, const container &fm, const container &fp, double beta, container &g)
Compute average along a fieldline .
Definition ds.h:981
void dss_centered(const FieldAligned &fa, double alpha, const container &fm, const container &f, const container &fp, double beta, container &g)
Centered derivative .
Definition ds.h:864
void ds_backward2(const FieldAligned &fa, double alpha, const container &fmm, const container &fm, const container &f, double beta, container &g)
2nd order backward derivative
Definition ds.h:813
void ds_divCentered(const FieldAligned &fa, double alpha, const container &fm, const container &fp, double beta, container &g)
centered derivative
Definition ds.h:951
void ds_forward2(const FieldAligned &fa, double alpha, const container &f, const container &fp, const container &fpp, double beta, container &g)
2nd order forward derivative
Definition ds.h:766
ONE FullLimiter
Full Limiter means there is a limiter everywhere.
Definition fieldaligned.h:30
void ds_forward(const FieldAligned &fa, double alpha, const container &f, const container &fp, double beta, container &g)
forward derivative
Definition ds.h:744
void ds_slope(const FieldAligned &fa, double alpha, const container &fm, const container &fp, double beta, container &g)
Compute simple slope along a fieldline .
Definition ds.h:1006
void swap_bc_perp(const FieldAligned &fa, const container &fm, const container &fp, container &fmg, container &fpg)
Swap the perp boundary condition.
Definition ds.h:307
@ einsPlus
plus interpolation in next plane
Definition fieldaligned.h:17
@ zeroForw
from dg to transformed coordinates
Definition fieldaligned.h:25
@ einsMinus
minus interpolation in previous plane
Definition fieldaligned.h:19
NotATensorTag tensor_category
This struct bundles a vector field and its divergence.
Definition fluxfunctions.h:444
Class for the evaluation of parallel derivatives.
Definition ds.h:349
void backward(double alpha, const container &f, double beta, container &g)
backward derivative
Definition ds.h:454
const FA & fieldaligned() const
Definition ds.h:685
void div(dg::direction dir, double alpha, const container &f, double beta, container &g)
Discretizes .
void backward(const container &f, container &g)
backward derivative
Definition ds.h:503
void divCentered(const container &f, container &g)
centered divergence
Definition ds.h:564
void divCentered(double alpha, const container &f, double beta, container &g)
centered divergence
Definition ds.h:544
void forward(const container &f, container &g)
forward derivative
Definition ds.h:512
void symv(double alpha, const container &f, double beta, container &g)
Discretizes as a symmetric matrix.
dg::geo::Fieldaligned< ProductGeometry, IMatrix, container > FA
Definition ds.h:350
void divForward(const container &f, container &g)
forward divergence
Definition ds.h:552
void divBackward(double alpha, const container &f, double beta, container &g)
backward divergence
Definition ds.h:536
DS(const dg::geo::CylindricalVectorLvl1 &vec, const ProductGeometry &grid, dg::bc bcx=dg::NEU, dg::bc bcy=dg::NEU, Limiter limit=FullLimiter(), double eps=1e-5, unsigned mx=10, unsigned my=10, double deltaPhi=-1, std::string interpolation_method="dg", bool benchmark=true)
Use the given vector field to construct.
Definition ds.h:382
void centered_bc_along_field(double alpha, const container &f, double beta, container &g, dg::bc bound, std::array< double, 2 > boundary_value={0, 0})
Same as dg::geo::ds_centered after dg::geo::ds_assign_bc_along_field_2nd.
Definition ds.h:487
void dss(const container &f, container &g)
Discretizes .
Definition ds.h:629
const container & weights() const
The volume form with dG weights.
Definition ds.h:675
DS(const dg::geo::TokamakMagneticField &vec, const ProductGeometry &grid, dg::bc bcx=dg::NEU, dg::bc bcy=dg::NEU, Limiter limit=FullLimiter(), double eps=1e-5, unsigned mx=10, unsigned my=10, double deltaPhi=-1, std::string interpolation_method="dg", bool benchmark=true)
Create the magnetic unit vector field and construct.
Definition ds.h:362
DS(FA fieldaligned)
Re-construct from a given Fieldaligned object.
void set_boundaries(dg::bc bcz, double left, double right)
Set boundary conditions in the limiter region.
Definition ds.h:413
void dssd(double alpha, const container &f, double beta, container &g)
Same as dg::geo::dssd_centered.
Definition ds.h:655
void construct(Params &&...ps)
Perfect forward parameters to one of the constructors.
Definition ds.h:406
void centered(double alpha, const container &f, double beta, container &g)
centered derivative
Definition ds.h:481
void set_boundaries(dg::bc bcz, const container &left, const container &right)
Set boundary conditions in the limiter region.
Definition ds.h:417
void dss(double alpha, const container &f, double beta, container &g)
Discretizes .
Definition ds.h:637
DS()
No memory allocation; all member calls except construct are invalid.
Definition ds.h:352
void dssd_bc_along_field(double alpha, const container &f, double beta, container &g, dg::bc bound, std::array< double, 2 > boundary_value={0, 0})
Same as dg::geo::dssd_centered after dg::geo::ds_assign_bc_along_field_2nd.
Definition ds.h:663
void symv(const container &f, container &g)
Discretizes .
Definition ds.h:612
void forward(double alpha, const container &f, double beta, container &g)
forward derivative
Definition ds.h:431
FA & fieldaligned()
access the underlying Fieldaligned object
Definition ds.h:684
void ds(dg::direction dir, double alpha, const container &f, double beta, container &g)
Discretizes .
void backward2(double alpha, const container &f, double beta, container &g)
2nd order backward derivative
Definition ds.h:465
void divForward(double alpha, const container &f, double beta, container &g)
forward divergence
Definition ds.h:528
void set_boundaries(dg::bc bcz, const container &global, double scal_left, double scal_right)
Set boundary conditions in the limiter region.
Definition ds.h:421
void centered(const container &f, container &g)
centered derivative
Definition ds.h:521
void div(dg::direction dir, const container &f, container &g)
Discretizes .
Definition ds.h:591
void forward2(double alpha, const container &f, double beta, container &g)
2nd order forward derivative
Definition ds.h:442
void dss_bc_along_field(double alpha, const container &f, double beta, container &g, dg::bc bound, std::array< double, 2 > boundary_value={0, 0})
Same as dg::geo::dss_centered after dg::geo::ds_assign_bc_along_field_2nd.
Definition ds.h:644
void ds(dg::direction dir, const container &f, container &g)
Discretizes .
Definition ds.h:574
void divBackward(const container &f, container &g)
backward divergence
Definition ds.h:558
Create and manage interpolation matrices from fieldline integration.
Definition fieldaligned.h:271
A tokamak field as given by R0, Psi and Ipol plus Meta-data like shape and equilibrium.
Definition magnetic_field.h:165