OSHUN  beta
Arbitrary Order Spherical-Harmonic 1D-3P Vlasov-Fokker-Planck-Maxwell code
Public Member Functions | Private Attributes
Magnetic_Field_1D_f1 Class Reference
Vlasov-Fokker-Planck operators in 1D

#include <vlasov_f1.h>

Collaboration diagram for Magnetic_Field_1D_f1:
Collaboration graph
[legend]

Public Member Functions

 Magnetic_Field_1D_f1 (size_t Nl, size_t Nm, double pmin, double pmax, size_t Np, double xmin, double xmax, size_t Nx)
 
void operator() (const DistFunc1D &Din, const Field1D &FBx, const Field1D &FBy, const Field1D &FBz, DistFunc1D &Dh)
 
void implicit (DistFunc1D &Din, const Field1D &FBx, const Field1D &FBy, const Field1D &FBz, double dt)
 

Private Attributes

SHarmonic1D FLM
 
valarray< complex< double > > A1
 
valarray< complex< double > > B1
 
Array2D< complex< double > > A2
 
complex< double > A3
 

Detailed Description

Definition at line 75 of file vlasov_f1.h.

Constructor & Destructor Documentation

◆ Magnetic_Field_1D_f1()

Magnetic_Field_1D_f1::Magnetic_Field_1D_f1 ( size_t  Nl,
size_t  Nm,
double  pmin,
double  pmax,
size_t  Np,
double  xmin,
double  xmax,
size_t  Nx 
)

Definition at line 298 of file vlasov_f1.cpp.

References A1, A3, and B1.

304  : A1(Nm+1), B1(Nl+1), A2(Nl+1,Nm+1), A3(0.5),
305  FLM(Np,Nx)
306 {
307 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
308  complex<double> lc, mc;
309  complex<double> c01(0.0,1.0);
310 
311 // Calculate the "A1" parameters
312 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
313  for (size_t m(0); m < Nm+1; ++m){
314  mc = static_cast< complex<double> >(m);
315  A1[m] = (-1.0)*c01*mc;
316  }
317 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
318 
319 // Calculate the "A2" parameters
320 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
321  // for (size_t l(1); l < Nl+1; ++l)
322  // for (size_t m=0; m<((Nm<l)?Nm:l)+1; ++m){
323  // lc = static_cast< complex<double> >(l);
324  // mc = static_cast< complex<double> >(m);
325  // A2(l,m) = (-0.5)*(lc+1.0-mc)*(lc+mc);
326  // }
327 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
328 
329 // Calculate the "A3" parameters
330 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
331  A3 = 0.5;
332 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
333 
334 // Calculate the "B1" parameters
335 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
336  for (size_t l(0); l < Nl+1; ++l){
337  lc = static_cast< complex<double> >(l);
338  B1[l] = (-1.0)*lc*(lc+1.0);
339  }
340 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
341 
342 }
Magnetic_Field_1D_f1::B1
valarray< complex< double > > B1
Definition: vlasov_f1.h:94
Magnetic_Field_1D_f1::FLM
SHarmonic1D FLM
Definition: vlasov_f1.h:92
Magnetic_Field_1D_f1::A1
valarray< complex< double > > A1
Definition: vlasov_f1.h:94
Magnetic_Field_1D_f1::A2
Array2D< complex< double > > A2
Definition: vlasov_f1.h:95
Magnetic_Field_1D_f1::A3
complex< double > A3
Definition: vlasov_f1.h:96

Member Function Documentation

◆ implicit()

void Magnetic_Field_1D_f1::implicit ( DistFunc1D &  Din,
const Field1D &  FBx,
const Field1D &  FBy,
const Field1D &  FBz,
double  dt 
)

Distribution function vector. 1 Nmx1 vector per momentum cell

Multiply Right Side to create right side vector

Unpack

Definition at line 422 of file vlasov_f1.cpp.

References Field1D::array(), DistFunc1D::q(), and Thomas_Tridiagonal().

424  {
425 //--------------------------------------------------------------
426 // This is the core calculation for the magnetic field
427 //--------------------------------------------------------------
428 
429  complex<double> ii(0.0,1.0);
430 
431  valarray<complex<double> > Bx(FBx.array());
432  valarray<complex<double> > Bm(FBy.array());
433  Bm *= (-1.0)*ii;
434  Bm += FBz.array();
435  valarray<complex<double> > Bp(FBy.array());
436  Bp *= ii;
437  Bp += FBz.array();
438 
439  Bx *= Din.q();//Din.mass();
440  Bm *= Din.q();//Din.mass();
441  Bp *= Din.q();//Din.mass();
442 
443  size_t l0(1);
444  size_t m0(1);
445 
446  int Nm(0);
447  size_t ind(0);
448  size_t indp(0);
449 
450  complex<double> temp(0.0,0.0);
451 
452  for (size_t i(0); i<Din(0,0).numx(); ++i)
453  {
454  for (size_t l(1); l < l0+1; ++l)
455  {
456  Nm = (l<m0?l:m0);
457  valarray<complex<double>> fin((0.0,0.0),2*Nm+1);
458  valarray<complex<double>> RHSv((0.0,0.0),2*Nm+1);
459  Array2D<complex<double>> LHS(2*Nm+1,2*Nm+1);
460  Array2D<complex<double>> RHS(2*Nm+1,2*Nm+1);
461 
462  //* The matrices corresponding to RHS and LHS are computed
463  // They apply to the vector
464  // [ f_l^m ] -> 0 index but applies to Nm th harmonic
465  // ...
466  // [ f_l^2 ]
467  // [ f_l^1 ]
468  // [ f_l^0 ] -> Nm index but applies to 0th harmonic
469  // [ f_l^-1 ]
470  // [ f_l^-2 ]
471  // ...
472  // [ f_l^-m ] *//
473  RHS(0,0) = 1.0-ii*double(Nm)*dt/2.0*Bx[i];
474  LHS(0,0) = conj(RHS(0,0));
475  RHS(0,1) = 0.25*dt*Bp[i];
476  LHS(0,1) = -0.25*dt*Bp[i];
477 
478  for (int m(Nm-1);m>0;--m)
479  {
480  ind = Nm-m;
481 
482  RHS(ind,ind) = 1.0-ii*double(m)*dt/2.0*Bx[i];
483  LHS(ind,ind) = conj(RHS(ind,ind));
484 
485  RHS(ind,ind+1) = 0.25*dt*Bp[i];
486  LHS(ind,ind+1) = -0.25*dt*Bp[i];
487 
488  RHS(ind,ind-1) = -0.25*dt*(l-m)*(l+m+1)*Bm[i];
489  LHS(ind,ind-1) = 0.25*dt*(l-m)*(l+m+1)*Bm[i];
490  }
491 
492  RHS(Nm,Nm) = 1.0;
493  LHS(Nm,Nm) = 1.0;
494  RHS(Nm,Nm-1) = -0.25*dt*(l)*(l+1)*Bm[i];
495  LHS(Nm,Nm-1) = 0.25*dt*(l)*(l+1)*Bm[i];
496 
497  //* Fill in bottom half which is a series of reflections and complex conjugates
498  RHS(Nm,Nm+1) = conj(RHS(Nm,Nm-1));
499  LHS(Nm,Nm+1) = conj(LHS(Nm,Nm-1));
500 
501  for (int m(-1);m>-Nm;--m)
502  {
503  ind = Nm-m;
504  indp = Nm-abs(m);
505 
506  RHS(ind,ind) = conj(RHS(indp,indp));
507  LHS(ind,ind) = conj(LHS(indp,indp));
508 
509  RHS(ind,ind+1) = conj(RHS(indp,indp-1));
510  LHS(ind,ind+1) = conj(LHS(indp,indp-1));
511 
512  RHS(ind,ind-1) = conj(RHS(indp,indp+1));
513  LHS(ind,ind-1) = conj(LHS(indp,indp+1));
514  }
515 
516  RHS(2*Nm,2*Nm) = conj(RHS(0,0));
517  LHS(2*Nm,2*Nm) = conj(LHS(0,0));
518  RHS(2*Nm,2*Nm-1) = conj(RHS(0,1));
519  LHS(2*Nm,2*Nm-1) = conj(LHS(0,1));
520 
521  for (size_t k(0); k < Din(0,0).nump(); ++k) {
523  fin[0] = Din(l, Nm)(k, i);
524  for (int m(Nm - 1); m > 1; --m) {
525  ind = Nm - m;
526  fin[ind] = Din(l, m)(k, i);
527  }
528  fin[Nm] = Din(l, 0)(k, i);
529  for (int m(-1); m > -Nm; --m) {
530  ind = Nm - m;
531  fin[ind] = conj(Din(l, abs(m))(k, i));
532  }
533  fin[2 * Nm] = conj(Din(l, Nm)(k, i));
534 
536  RHSv[0] = fin[0] * RHS(0, 0) + fin[1] * RHS(0, 1);
537  for (size_t mm(1); mm < 2 * Nm; ++mm) {
538  RHSv[mm] = fin[mm - 1] * RHS(mm, mm - 1) + fin[mm] * RHS(mm, mm) + fin[mm + 1] * RHS(mm, mm + 1);
539  }
540  RHSv[2 * Nm] = fin[2 * Nm - 1] * RHS(2 * Nm, 2 * Nm - 1) + fin[2 * Nm] * RHS(2 * Nm, 2 * Nm);
541 
542 // std::cout << "\n\n LHS = \n";
543 // for (size_t i(0); i < LHS.dim1(); ++i) {
544 // std::cout << "i = " << i << " :::: ";
545 // for (size_t j(0); j < LHS.dim2(); ++j) {
546 // std::cout << LHS(i, j) << " ";
547 // }
548 // std::cout << "\n";
549 // }
550 
551  Thomas_Tridiagonal(LHS,RHSv,fin);
552 
554  Din(l,Nm)(k,i) = fin[0];
555  for (int m(Nm-1);m>1;--m)
556  {
557  ind = Nm-m;
558  Din(l,m)(k,i) = fin[ind];
559  }
560  Din(l,0)(k,i) = fin[Nm];
561 
562  }
563 
564  }
565  }
566 
567 }
Thomas_Tridiagonal
bool Thomas_Tridiagonal(Array2D< double > &A, valarray< double > &d, valarray< double > &xk)
The tridiagonal solver for implicit collisions.
Definition: nmethods.cpp:216
Field1D::array
valarray< complex< double > > & array() const
Definition: state.h:196
DistFunc1D::q
double q() const
Definition: state.h:399
Here is the call graph for this function:

◆ operator()()

void Magnetic_Field_1D_f1::operator() ( const DistFunc1D &  Din,
const Field1D &  FBx,
const Field1D &  FBy,
const Field1D &  FBz,
DistFunc1D &  Dh 
)

Definition at line 347 of file vlasov_f1.cpp.

References A1, A3, Field1D::array(), B1, FLM, SHarmonic1D::mxaxis(), DistFunc1D::q(), and SHarmonic1D::Re().

349  {
350 //--------------------------------------------------------------
351 // This is the core calculation for the magnetic field
352 //--------------------------------------------------------------
353 
354  complex<double> ii(0.0,1.0);
355 
356  valarray<complex<double> > Bx(FBx.array());
357  valarray<complex<double> > Bm(FBy.array());
358  Bm *= (-1.0)*ii;
359  Bm += FBz.array();
360  valarray<complex<double> > Bp(FBy.array());
361  Bp *= ii;
362  Bp += FBz.array();
363 
364  Bx *= Din.q();
365  Bm *= Din.q();
366  Bp *= Din.q();
367 
368  size_t l0(B1.size()-1);
369  size_t m0(A1.size()-1);
370 
371 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
372 // m = 0, 1 < l < l0+1
373 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
374  Bp *= A3;
375  for (size_t l(1); l < l0+1; ++l){
376  FLM = Din(l,0); Dh(l,1) += FLM.mxaxis(Bp);
377  }
378 
379 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
380 // m = 1, l = 1
381 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
382  FLM = Din(1,1); Bx *= A1[1]; Dh(1,1) += FLM.mxaxis(Bx);
383  FLM = Din(1,1); Bm *= B1[1]; FLM = FLM.mxaxis(Bm); Dh(1,0) += FLM.Re();
384 
385 // // - - - - - - - - - - - - - - - - - - - - - - - - - - -
386 // // m = 1, l > 1
387 // // - - - - - - - - - - - - - - - - - - - - - - - - - - -
388 // for (size_t l(2); l < l0+1; ++l){
389 // FLM = Din(l,1); Dh(l,2) += FLM.mxaxis(Bp);
390 // FLM = Din(l,1); Dh(l,1) += FLM.mxaxis(Bx);
391 // FLM = Din(l,1); Bm *= B1[l]/B1[l-1]; FLM = FLM.mxaxis(Bm); Dh(l,0) += FLM.Re();
392 // }
393 // Bm *= 1.0/B1[l0];
394 
395 // // - - - - - - - - - - - - - - - - - - - - - - - - - - -
396 // // m > 1, l = m
397 // // - - - - - - - - - - - - - - - - - - - - - - - - - - -
398 // for (size_t m(2); m < m0; ++m){
399 // FLM = Din(m,m); Bx *= A1[m]/A1[m-1]; Dh(m,m ) += FLM.mxaxis(Bx);
400 // FLM = Din(m,m); Bm *= A2(m,m); Dh(m,m-1) += FLM.mxaxis(Bm);
401 // for (size_t l(m+1); l < l0+1; ++l){
402 // FLM = Din(l,m); Dh(l,m+1) += FLM.mxaxis(Bp);
403 // FLM = Din(l,m); Dh(l,m ) += FLM.mxaxis(Bx);
404 // FLM = Din(l,m); Bm *= A2(l,m)/A2(l-1,m); Dh(l,m-1) += FLM.mxaxis(Bm);
405 // }
406 // Bm *= 1.0/A2(l0,m);
407 // }
408 
409 // // - - - - - - - - - - - - - - - - - - - - - - - - - - -
410 // // m = m0, l >= m0
411 // // - - - - - - - - - - - - - - - - - - - - - - - - - - -
412 // FLM = Din(m0,m0); Bx *= A1[m0]/A1[m0-1]; Dh(m0,m0) += FLM.mxaxis(Bx);
413 // FLM = Din(m0,m0); Bm *= A2(m0,m0)/*/A2(l0,m0-1)*/; Dh(m0,m0-1) += FLM.mxaxis(Bm);
414 // for (size_t l(m0+1); l < l0+1; ++l){
415 // FLM = Din(l,m0); Dh(l,m0 ) += FLM.mxaxis(Bx);
416 // FLM = Din(l,m0); Bm *= A2(l,m0)/A2(l-1,m0); Dh(l,m0-1) += FLM.mxaxis(Bm);
417 // }
418 
419 }
Magnetic_Field_1D_f1::B1
valarray< complex< double > > B1
Definition: vlasov_f1.h:94
Magnetic_Field_1D_f1::FLM
SHarmonic1D FLM
Definition: vlasov_f1.h:92
Magnetic_Field_1D_f1::A1
valarray< complex< double > > A1
Definition: vlasov_f1.h:94
Magnetic_Field_1D_f1::A3
complex< double > A3
Definition: vlasov_f1.h:96
SHarmonic1D::mxaxis
SHarmonic1D & mxaxis(const valarray< complex< double > > &shmulti)
Definition: state.cpp:126
SHarmonic1D::Re
SHarmonic1D & Re()
Definition: state.cpp:130
Field1D::array
valarray< complex< double > > & array() const
Definition: state.h:196
DistFunc1D::q
double q() const
Definition: state.h:399
Here is the call graph for this function:

Field Documentation

◆ A1

valarray< complex<double> > Magnetic_Field_1D_f1::A1
private

Definition at line 94 of file vlasov_f1.h.

Referenced by Magnetic_Field_1D_f1(), and operator()().

◆ A2

Array2D< complex<double> > Magnetic_Field_1D_f1::A2
private

Definition at line 95 of file vlasov_f1.h.

◆ A3

complex<double> Magnetic_Field_1D_f1::A3
private

Definition at line 96 of file vlasov_f1.h.

Referenced by Magnetic_Field_1D_f1(), and operator()().

◆ B1

valarray< complex<double> > Magnetic_Field_1D_f1::B1
private

Definition at line 94 of file vlasov_f1.h.

Referenced by Magnetic_Field_1D_f1(), and operator()().

◆ FLM

SHarmonic1D Magnetic_Field_1D_f1::FLM
private

Definition at line 92 of file vlasov_f1.h.

Referenced by operator()().

The documentation for this class was generated from the following files: