Extension: Geometries
#include "dg/geometries/geometries.h"
Loading...
Searching...
No Matches
ds.h
Go to the documentation of this file.
1#pragma once
2
3#include "dg/algorithm.h"
4#include "fieldaligned.h"
5#ifdef MPI_VERSION
6#include "mpi_fieldaligned.h"
7#endif //MPI_VERSION
8#include "magnetic_field.h"
9
15namespace dg{
16namespace geo{
17
65namespace detail
66{
67struct DSCentered
68{
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,
72 double hp)
73 {
74 dsf = m_alpha*(
75 fm*( 1./(hp+hm) - 1./hm) +
76 fo*( 1./hm - 1./hp) +
77 fp*( 1./hp - 1./(hp+hm))
78 ) + m_beta*dsf;
79 };
80
81 private:
82 double m_alpha, m_beta;
83};
84struct DSSCentered
85{
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,
89 double hp)
90 {
91 dssf = m_alpha*(
92 2.*fm/(hp+hm)/hm - 2.*fo/hp/hm + 2.*fp/(hp+hm)/hp
93 ) + m_beta*dssf;
94 };
95
97 void operator()( double& dssf, double fm, double fo, double fp,
98 double bPm, double bP0, double bPp)
99 {
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;
104
105 dssf = m_alpha*bP0*( bP2*fp2 - bM2*fm2)/m_delta + m_beta*dssf;
106 }
107
108 private:
109 double m_alpha, m_beta, m_delta;
110};
111struct DSSDCentered
112{
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)
118 {
119 // various combinations of bP do not seem to matter
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.;
126
127 dssdf = m_alpha*( gp2*fp2*bP2*bP2 - bM2*bM2*gm2*fm2)/m_delta + m_beta*dssdf;
128
129 // does not seem to conserve nicely
130 //dssdf = m_alpha*( fp*gp2*bP2*bP2 - 2*bP0*bP0*fo + bM2*bM2*gm2*fm)/(hp)/(hm) + m_beta*dssdf;
131 //dssdf = m_alpha*( fp*Gp/Go*bPp*bPp - 2*bP0*bP0*fo + bPm*bPm*Gm/Go*fm)/(hp)/(hm) + m_beta*dssdf;
132 };
133
134 private:
135 double m_alpha, m_beta, m_delta;
136};
137
138}//namespace detail
140
169template<class FieldAligned, class container>
170void assign_bc_along_field_2nd( const FieldAligned& fa, const container& fm,
171 const container& f, const container& fp, container& fmg, container& fpg,
172 dg::bc bound, std::array<double,2> boundary_value = {0,0})
173{
174 double delta = fa.deltaPhi();
175 if( bound == dg::NEU)
176 {
177 double dbm = boundary_value[0], dbp = boundary_value[1];
178 dg::blas1::subroutine( [dbm, dbp, delta]DG_DEVICE( double fm, double fo,
179 double fp, double& fmg, double& fpg,
180 double hbm, double hbp, double bbm, double bbo, double bbp
181 ){
182 //formula derived for non-equidistant grid
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() );
199 }
200 else// if( bound == dg:DIR)
201 {
202 double fbm = boundary_value[0], fbp = boundary_value[1];
203 dg::blas1::subroutine( [fbm, fbp, delta]DG_DEVICE( double fm, double fo,
204 double fp, double& fmg, double& fpg,
205 double hbm, double hbp, double bbm, double bbo, double bbp
206 ){
207 //formula derived for non-equidistant grid
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(),
221 fa.bbo(), fa.bbp());
222 }
223
224}
241template<class FieldAligned, class container>
242void assign_bc_along_field_1st( const FieldAligned& fa, const container& fm,
243 const container& fp, container& fmg, container& fpg,
244 dg::bc bound, std::array<double,2> boundary_value = {0,0})
245{
246 double delta = fa.deltaPhi();
247 if( bound == dg::NEU)
248 {
249 double dbm = boundary_value[0], dbp = boundary_value[1];
250 dg::blas1::subroutine( [dbm, dbp, delta]DG_DEVICE( double fm, double fp,
251 double& fmg, double& fpg, double bbm, double bbp
252 ){
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() );
260 }
261 else// if( bound == dg:DIR)
262 {
263 double fbm = boundary_value[0], fbp = boundary_value[1];
264 dg::blas1::subroutine( [fbm, fbp, delta]DG_DEVICE( double fm, double fp,
265 double& fmg, double& fpg, double hbm,
266 double hbp, double bbm, double bbo, double bbp
267 ){
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(),
277 fa.bbo(), fa.bbp());
278 }
279}
280
306template<class FieldAligned, class container>
307void swap_bc_perp( const FieldAligned& fa, const container& fm,
308 const container& fp, container& fmg, container& fpg)
309{
310 dg::blas1::subroutine( []DG_DEVICE( double fm, double fp,
311 double& fmg, double& fpg,
312 double bbm, double bbo, double bbp
313 ){
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() );
317
318}
319
347template< class ProductGeometry, class IMatrix, class container >
348struct DS
349{
352 DS(){}
353
361 template <class Limiter>
362 DS(const dg::geo::TokamakMagneticField& vec, const ProductGeometry& grid,
363 dg::bc bcx = dg::NEU,
364 dg::bc bcy = dg::NEU,
365 Limiter limit = FullLimiter(),
366 double eps = 1e-5,
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) )
372 {
373 }
381 template<class Limiter>
382 DS(const dg::geo::CylindricalVectorLvl1& vec, const ProductGeometry& grid,
383 dg::bc bcx = dg::NEU,
384 dg::bc bcy = dg::NEU,
385 Limiter limit = FullLimiter(),
386 double eps = 1e-5,
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))
392 {
393 }
405 template<class ...Params>
406 void construct( Params&& ...ps)
407 {
408 //construct and swap
409 *this = DS( std::forward<Params>( ps)...);
410 }
411
413 void set_boundaries( dg::bc bcz, double left, double right){
414 m_fa.set_boundaries( bcz, left, right);
415 }
417 void set_boundaries( dg::bc bcz, const container& left, const container& right){
418 m_fa.set_boundaries( bcz, left, right);
419 }
421 void set_boundaries( dg::bc bcz, const container& global, double scal_left, double scal_right){
422 m_fa.set_boundaries( bcz, global, scal_left, scal_right);
423 }
424
431 void forward( double alpha, const container& f, double beta, container& g){
432 m_fa(zeroForw, f, m_tempO);
433 m_fa(einsPlus, f, m_tempP);
434 ds_forward( m_fa, alpha, m_tempO, m_tempP, beta, g);
435 }
442 void forward2( double alpha, const container& f, double beta, container& g){
443 m_fa(zeroForw, f, m_tempO);
444 m_fa(einsPlus, f, m_tempP);
445 m_fa(einsPlus, m_tempP, m_tempM);
446 ds_forward2( m_fa, alpha, m_tempO, m_tempP, m_tempM, beta, g);
447 }
454 void backward( double alpha, const container& f, double beta, container& g){
455 m_fa(einsMinus, f, m_tempM);
456 m_fa(zeroForw, f, m_tempO);
457 ds_backward( m_fa, alpha, m_tempM, m_tempO, beta, g);
458 }
465 void backward2( double alpha, const container& f, double beta, container& g){
466 m_fa(einsMinus, f, m_tempM);
467 m_fa(einsMinus, m_tempM, m_tempP);
468 m_fa(zeroForw, f, m_tempO);
469 ds_backward2( m_fa, alpha, m_tempP, m_tempM, m_tempO, beta, g);
470 }
481 void centered( double alpha, const container& f, double beta, container& g){
482 m_fa(einsPlus, f, m_tempP);
483 m_fa(einsMinus, f, m_tempM);
484 ds_centered( m_fa, alpha, m_tempM, m_tempP, beta, g);
485 }
488 double alpha, const container& f, double beta, container& g, dg::bc bound,
489 std::array<double,2> boundary_value = {0,0}){
490 m_fa(einsPlus, f, m_tempP);
491 m_fa(einsMinus, f, m_tempM);
492 assign_bc_along_field_2nd( m_fa, m_tempM, f, m_tempP, m_tempM, m_tempP,
493 bound, boundary_value);
494 ds_centered( m_fa, alpha, m_tempM, m_tempP, beta, g);
495
496 }
503 void backward( const container& f, container& g){
504 backward(1., f,0.,g);
505 }
512 void forward( const container& f, container& g){
513 forward(1.,f, 0.,g);
514 }
521 void centered( const container& f, container& g){
522 centered(1.,f,0.,g);
523 }
524
528 void divForward( double alpha, const container& f, double beta, container& g){
529 m_fa(einsPlus, f, m_tempP);
530 m_fa(zeroForw, f, m_tempO);
531 ds_divForward( m_fa, alpha, m_tempO, m_tempP, beta, g);
532 }
536 void divBackward( double alpha, const container& f, double beta, container& g){
537 m_fa(einsMinus, f, m_tempM);
538 m_fa(zeroForw, f, m_tempO);
539 ds_divBackward( m_fa, alpha, m_tempM, m_tempO, beta, g);
540 }
544 void divCentered(double alpha, const container& f, double beta, container& g){
545 m_fa(einsPlus, f, m_tempP);
546 m_fa(einsMinus, f, m_tempM);
547 ds_divCentered( m_fa, alpha, m_tempM, m_tempP, beta, g);
548 }
552 void divForward(const container& f, container& g){
553 divForward( 1.,f,0.,g);
554 }
558 void divBackward(const container& f, container& g){
559 divBackward( 1.,f,0.,g);
560 }
564 void divCentered(const container& f, container& g){
565 divCentered( 1.,f,0.,g);
566 }
567
574 void ds(dg::direction dir, const container& f, container& g){
575 ds(dir, 1., f, 0., g);
576 }
583 void ds(dg::direction dir, double alpha, const container& f, double beta, container& g);
591 void div(dg::direction dir, const container& f, container& g){
592 div(dir, 1., f, 0., g);
593 }
601 void div(dg::direction dir, double alpha, const container& f, double beta, container& g);
602
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){
630 dss( 1., f, 0., g);}
637 void dss( double alpha, const container& f, double beta, container& g){
638 m_fa(einsPlus, f, m_tempP);
639 m_fa(einsMinus, f, m_tempM);
640 m_fa(zeroForw, f, m_tempO);
641 dss_centered( m_fa, alpha, m_tempM, m_tempO, m_tempP, beta, g);
642 }
645 double alpha, const container& f, double beta, container& g, dg::bc bound,
646 std::array<double,2> boundary_value = {0,0}){
647 m_fa(einsPlus, f, m_tempP);
648 m_fa(einsMinus, f, m_tempM);
649 m_fa(zeroForw, f, m_tempO);
650 assign_bc_along_field_2nd( m_fa, m_tempM, m_tempO, m_tempP, m_tempM, m_tempP,
651 bound, boundary_value);
652 dss_centered( m_fa, alpha, m_tempM, m_tempO, m_tempP, beta, g);
653 }
655 void dssd( double alpha, const container& f, double
656 beta, container& g){
657 m_fa(einsPlus, f, m_tempP);
658 m_fa(einsMinus, f, m_tempM);
659 m_fa(zeroForw, f, m_tempO);
660 dssd_centered( m_fa, alpha, m_tempM, m_tempO, m_tempP, beta, g);
661 }
663 void dssd_bc_along_field( double alpha, const
664 container& f, double beta, container& g, dg::bc bound,
665 std::array<double,2> boundary_value = {0,0}){
666 m_fa(einsPlus, f, m_tempP);
667 m_fa(einsMinus, f, m_tempM);
668 m_fa(zeroForw, f, m_tempO);
669 assign_bc_along_field_2nd( m_fa, m_tempM, m_tempO, m_tempP, m_tempM, m_tempP,
670 bound, boundary_value);
671 dssd_centered( m_fa, alpha, m_tempM, f, m_tempP, beta, g);
672 }
673
675 const container& weights()const {
676 return m_fa.sqrtG();
677 }
678
684 FA& fieldaligned(){return m_fa;}
685 const FA& fieldaligned()const{return m_fa;}
686 private:
688 container m_tempP, m_tempO, m_tempM;
689};
690
693
694template<class Geometry, class I, class container>
695DS<Geometry, I, container>::DS( Fieldaligned<Geometry, I, container> fa): m_fa(fa)
696{
697 m_tempP = fa.sqrtG(), m_tempM = m_tempO = m_tempP;
698}
699
700template<class G, class I, class container>
701inline void DS<G,I,container>::ds( dg::direction dir, double alpha,
702 const container& f, double beta, container& dsf) {
703 switch( dir){
704 case dg::centered:
705 return centered( alpha, f, beta, dsf);
706 case dg::forward:
707 return forward( alpha, f, beta, dsf);
708 case dg::backward:
709 return backward( alpha, f, beta, dsf);
710 }
711}
712template<class G, class I, class container>
713inline void DS<G,I,container>::div( dg::direction dir, double alpha,
714 const container& f, double beta, container& dsf) {
715 switch( dir){
716 case dg::centered:
717 return divCentered( alpha, f, beta, dsf);
718 case dg::forward:
719 return divForward( alpha, f, beta, dsf);
720 case dg::backward:
721 return divBackward( alpha, f, beta, dsf);
722 }
723}
724
725
726template<class G,class I, class container>
727void DS<G,I,container>::symv( double alpha, const container& f, double beta, container& dsTdsf)
728{
729 dssd( alpha, f, beta, dsTdsf);
730};
732
743template<class FieldAligned, class container>
744void ds_forward(const FieldAligned& fa, double alpha, const container& f,
745 const container& fp, double beta, container& g)
746{
747 //direct
748 double delta = fa.deltaPhi();
749 dg::blas1::subroutine( [ alpha, beta, delta]DG_DEVICE(
750 double& dsf, double fo, double fp, double bphi){
751 dsf = alpha*bphi*( fp - fo)/delta + beta*dsf;
752 },
753 g, f, fp, fa.bphi());
754}
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)
768{
769 //direct
770 double delta = fa.deltaPhi();
771 dg::blas1::subroutine( [ alpha, beta, delta]DG_DEVICE(
772 double& dsf, double fo, double fp, double fpp, double bphi){
773 dsf = alpha*bphi*( -3.*fo + 4.*fp - fpp)/2./delta
774 + beta*dsf;
775 },
776 g, f, fp, fpp, fa.bphi());
777}
778
789template<class FieldAligned, class container>
790void ds_backward( const FieldAligned& fa, double alpha, const container& fm,
791 const container& f, double beta, container& g)
792{
793 //direct
794 double delta = fa.deltaPhi();
795 dg::blas1::subroutine( [ alpha, beta, delta] DG_DEVICE(
796 double& dsf, double fo, double fm, double bphi){
797 dsf = alpha*bphi*( fo - fm)/delta + beta*dsf;
798 },
799 g, f, fm, fa.bphi());
800
801}
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)
814{
815 //direct
816 double delta = fa.deltaPhi();
817 dg::blas1::subroutine( [ alpha, beta, delta] DG_DEVICE(
818 double& dsf, double fo, double fm, double fmm, double bphi){
819 dsf = alpha*bphi*( 3.*fo - 4.*fm + fmm)/2./delta
820 + beta*dsf;
821 },
822 g, f, fm, fmm, fa.bphi());
823
824}
825
826
841template<class FieldAligned, class container>
842void ds_centered( const FieldAligned& fa, double alpha, const container& fm,
843 const container& fp, double beta, container& g)
844{
845 //direct discretisation
846 double delta=fa.deltaPhi();
847 dg::blas1::subroutine( [alpha,beta,delta]DG_DEVICE( double& g, double fm,
848 double fp, double bphi){
849 g = alpha*bphi*(fp-fm)/2./delta + beta*g;
850 }, g, fm, fp, fa.bphi());
851}
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)
866{
867 dg::blas1::subroutine( detail::DSSCentered( alpha, beta, fa.deltaPhi()),
868 g, fm, f, fp, fa.bphiM(), fa.bphi(), fa.bphiP());
869}
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)
884{
885 dg::blas1::subroutine( detail::DSSDCentered( alpha, beta, fa.deltaPhi()),
886 g, fm, f, fp, fa.sqrtGm(), fa.sqrtG(), fa.sqrtGp(),
887 fa.bphiM(), fa.bphi(), fa.bphiP());
888}
889
903template<class FieldAligned, class container>
904void ds_divBackward( const FieldAligned& fa, double alpha, const container& fm,
905 const container& f, double beta, container& g)
906{
907 double delta = fa.deltaPhi();
908 dg::blas1::subroutine( [alpha,beta,delta] DG_DEVICE( double& dsf, double f0,
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());
912}
913
927template<class FieldAligned, class container>
928void ds_divForward( const FieldAligned& fa, double alpha, const container& f,
929 const container& fp, double beta, container& g)
930{
931 double delta = fa.deltaPhi();
932 dg::blas1::subroutine( [alpha,beta,delta] DG_DEVICE( double& dsf, double f0,
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());
936}
950template<class FieldAligned, class container>
951void ds_divCentered( const FieldAligned& fa, double alpha, const container& fm, const container& fp,
952 double beta, container& g)
953{
954 double delta = fa.deltaPhi();
955 dg::blas1::subroutine( [alpha,beta,delta]DG_DEVICE( double& dsf, double fm,
956 double fp, double Gm, double Gp, double G0,
957 double bPm, double bP0, double bPp)
958 {
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());
962
963}
964
980template<class FieldAligned, class container>
981void ds_average( const FieldAligned& fa, double alpha,
982 const container& fm, const container& fp, double beta, container& g)
983{
984 dg::blas1::subroutine( [alpha,beta]DG_DEVICE( double& g, double fm, double fp
985 ){
986 g = alpha*(fp+fm)/2. + beta*g;
987 }, g, fm, fp);
988}
1005template<class FieldAligned, class container>
1006void ds_slope( const FieldAligned& fa, double alpha,
1007 const container& fm, const container& fp, double beta, container& g)
1008{
1009 ds_centered( fa, alpha, fm, fp, beta, g);
1010}
1011
1012
1013}//namespace geo
1014
1016template< class G, class I, class V>
1017struct TensorTraits< geo::DS<G,I, V> >
1018{
1019 using value_type = double;
1020 using tensor_category = SelfMadeMatrixTag;
1021};
1023}//namespace dg
void plus(ContainerType &x, value_type alpha)
void subroutine(Subroutine f, ContainerType &&x, ContainerTypes &&... xs)
direction
backward
centered
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
#define DG_DEVICE
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