OSHUN  beta
Arbitrary Order Spherical-Harmonic 1D-3P Vlasov-Fokker-Planck-Maxwell code
vlasov_f1.cpp
Go to the documentation of this file.
1 
9 //--------------------------------------------------------------
10 // Standard libraries
11 #include <iostream>
12 #include <vector>
13 #include <valarray>
14 #include <complex>
15 #include <algorithm>
16 #include <cstdlib>
17 
18 #include <math.h>
19 #include <map>
20 
21 // My libraries
22 #include "lib-array.h"
23 #include "lib-algorithms.h"
24 #include "nmethods.h"
25 
26 // Declerations
27 #include "state.h"
28 #include "formulary.h"
29 #include "input.h"
30 #include "fluid.h"
31 #include "vlasov.h"
32 #include "vlasov_f1.h"
33 
34 
35 //**************************************************************
36 //--------------------------------------------------------------
37 Electric_Field_1D_f1::Electric_Field_1D_f1(size_t Nl, size_t Nm,
38  double pmin, double pmax, size_t Np,
39  double xmin, double xmax, size_t Nx)
40 //--------------------------------------------------------------
41 // Constructor
42 //--------------------------------------------------------------
43  : Hp0(Nl+1),
44  H(Np,Nx), G(Np,Nx), TMP(Np,Nx),
45  pr(Algorithms::MakeAxis(static_cast<complex<double> >(pmin),
46  static_cast<complex<double> >(pmax),
47  Np)),
48  invpr(pr)
49 {
50 
51 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
52  complex<double> lc, mc;
53 
54 // Inverted momentum axis
55  for (size_t i(0); i < pr.size(); ++i) {
56  invpr[i] = 1.0/pr[i];
57  }
58  double idp = (-1.0)/ (2.0*(pmax-pmin)/double(Np-1)); // -1/(2dp)
59 
60 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
61  A100 = static_cast<complex<double>>(idp);
62  C100 = static_cast<complex<double>>(0.5*idp);
63  A210 = static_cast<complex<double>>(1.0/3.0*idp);
64  B211 = static_cast<complex<double>>(2.0/3.0);
65  A310 = static_cast<complex<double>>(2.0/5.0*idp);
66  C311 = static_cast<complex<double>>(-0.5*idp/5.0);
67 
68 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
69 // -2*Dp*H at the 0 momentum cell
70  Hp0[0] = 1.0 / pr[0];
71  for (size_t l(1); l < Nl+1; ++l) {
72  double ld(l);
73  Hp0[l] = Hp0[l-1] * (pr[0]/pr[1]) * (2.0*ld+1.0)/(2.0*ld-1.0);
74  }
75  Hp0 *= (-2.0)*(pr[1]-pr[0]);
76 
77 }
78 //--------------------------------------------------------------
79 
80 //--------------------------------------------------------------
81 void Electric_Field_1D_f1::operator()(const DistFunc1D& Din,
82  const Field1D& FEx, const Field1D& FEy, const Field1D& FEz,
83  DistFunc1D& Dh) {
84 //--------------------------------------------------------------
85 // This is the core calculation for the electric field
86 //--------------------------------------------------------------
87 
88  complex<double> ii(0.0,1.0);
89 
90  valarray<complex<double> > Ex(FEx.array());
91  valarray<complex<double> > Em(FEz.array());
92  Em *= (-1.0)*ii;
93  Em += FEy.array();
94  valarray<complex<double> > Ep(FEz.array());
95  Ep *= ii;
96  Ep += FEy.array();
97 
98  Ex *= Din.q();
99  Em *= Din.q();
100  Ep *= Din.q();
101 
102  size_t l0(Din.l0());
103  size_t m0(Din.m0());
104 
105 
106 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107 // m = 0, l = 0
108 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
109  MakeG00(Din(0,0));
110  Ex *= A100; TMP = G; Dh(1,0) += G.mxaxis(Ex);
111  Em *= C100; Dh(1,1) += TMP.mxaxis(Em);
112 
113 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
114 // m = 0, l = 1
115 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
116  MakeGH(Din(1,0),1);
117  Ex *= A210 / A100; Dh(0,0) += H.mxaxis(Ex);
118 
119 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
120 // m = 0, 1 < l < l0
121 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
122 
123  // Din(2,0) = Din(0,0);
124  // Din(2,0) *= static_cast<complex<double>>(1.0/3.0);
125 
126  // MakeGH(Din(2,0),2);
127  // Ex *= A310 / A210; TMP = H; Dh(1,0) += H.mxaxis(Ex);
128  // Em *= C311 / C100; Dh(1,1) += TMP.mxaxis(Em);
129 
130 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
131 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
132 // m = 1, l = 1
133 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
134  MakeGH(Din(1,1),1);
135  Ep *= B211; H = H.mxaxis(Ep); Dh(0,0) += H.Re();
136 
137 }
138 //--------------------------------------------------------------
139 
140 
141 //--------------------------------------------------------------
142 void Electric_Field_1D_f1::Implicit_Ex(const DistFunc1D& Din, const Field1D& FEx, DistFunc1D& Dh) {
143 //--------------------------------------------------------------
144 // This is the calculation for the implicit electric field in x,
145 // which is to say Ex as it acts on the first few harmonics.
146 //--------------------------------------------------------------
147 
148  valarray<complex<double> > Ex(FEx.array());
149  Ex *= Din.q();
150 
151 // m = 0, l = 0
152  MakeG00(Din(0,0));
153  Ex *= A100; Dh(1,0) += G.mxaxis(Ex);
154 
155 // m = 0, l = 1
156  MakeGH(Din(1,0),1);
157  Ex *= A210 / A100; Dh(0,0) += H.mxaxis(Ex);
158  // Ex *= A1(1,0) / A2(1,0); Dh(2,0) += G.mxaxis(Ex);
159 
160 // m = 0, l = 2
161  // Din(2,0) = Din(0,0);
162  // Din(2,0) *= static_cast<complex<double>>(1.0/3.0);
163 
164  // MakeGH(Din(2,0),2);
165  // Ex *= A310 / A210; TMP = H; Dh(1,0) += H.mxaxis(Ex);
166 
167 
168 }
169 //--------------------------------------------------------------
170 //--------------------------------------------------------------
171 void Electric_Field_1D_f1::Implicit_Ey(const DistFunc1D& Din, const Field1D& FEy, DistFunc1D& Dh) {
172 //--------------------------------------------------------------
173 // This is the calculation for the implicit electric field in y
174 // and z, which is to say Ex as it acts on the first few harmonics.
175 //--------------------------------------------------------------
176 
177  valarray<complex<double> > Em(FEy.array());
178  valarray<complex<double> > Ep(FEy.array());
179 
180 
181  Em *= Din.q();;
182  Ep *= Din.q();;
183 
184 // m = 0, l = 0
185  MakeG00(Din(0,0));
186  Em *= C100; Dh(1,1) += G.mxaxis(Em);
187 
188 // m = 0, l = 1
189  // MakeGH(Din(1,0),1);
190  // Em *= C1[1] / C1[0]; Dh(2,1) += G.mxaxis(Em);
191 
192 // m = 0, 1 < l < l0
193  // Din(2,0) = Din(0,0);
194  // Din(2,0) *= static_cast<complex<double>>(1.0/3.0);
195 
196  // MakeGH(Din(2,0),2);
197  // Em *= C311 / C100; Dh(1,1) += TMP.mxaxis(Em);
198 
199 // m = 0, l = 3
200  // MakeGH(Din(3,0),3);
201  // Em *= C3[3] / C3[2]; Dh(2,1) += H.mxaxis(Em);
202 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
203 
204 
205 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
206 // m = 1, l = 1
207  MakeGH(Din(1,1),1);
208  Ep *= B211; H = H.mxaxis(Ep); Dh(0,0) += H.Re();
209 }
210 //
212 void Electric_Field_1D_f1::Implicit_Ez(const DistFunc1D& Din, const Field1D& FEz, DistFunc1D& Dh) {
213 //--------------------------------------------------------------
214 // This is the calculation for the implicit electric field in y
215 // and z, which is to say Ex as it acts on the first few harmonics.
216 //--------------------------------------------------------------
217  complex<double> ii(0.0,1.0);
218 
219 
220  valarray<complex<double> > Em(FEz.array());
221  Em *= (-1.0)*ii;
222 
223 
224  valarray<complex<double> > Ep(FEz.array());
225  Ep *= ii;
226  // m = 0, l = 0
227  MakeG00(Din(0,0));
228  Em *= C100; Dh(1,1) += G.mxaxis(Em);
229 
230 // m = 0, l = 1
231  // MakeGH(Din(1,0),1);
232  // Em *= C1[1] / C1[0]; Dh(2,1) += G.mxaxis(Em);
233 
234 // m = 0, 1 < l < l0
235  Din(2,0) = Din(0,0);
236  Din(2,0) *= static_cast<complex<double>>(1.0/3.0);
237 
238  MakeGH(Din(2,0),2);
239  Em *= C311 / C100; Dh(1,1) += TMP.mxaxis(Em);
240 
241 // // m = 0, l = 3
242 // MakeGH(Din(3,0),3);
243 // Em *= C3[3] / C3[2]; Dh(2,1) += H.mxaxis(Em);
244 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
245 
246 
247 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
248 // m = 1, l = 1
249  MakeGH(Din(1,1),1);
250  Ep *= B211; H = H.mxaxis(Ep); Dh(0,0) += H.Re();
251 }
252 //--------------------------------------------------------------
253 
254 
255 //--------------------------------------------------------------
256 // Make derivatives 2*Dp*(l+1/l)*G and -2*Dp*H for a given f
257 void Electric_Field_1D_f1::MakeGH(SHarmonic1D& f, size_t el){
258 //--------------------------------------------------------------
259  valarray<complex<double> > invpax(invpr);
260  complex<double> ld(el);
261 
262  invpax *= (-2.0)*(ld+1.0) * (pr[1]-pr[0]);
263 
264  G = f; H = f;
265  G = G.Dp();
266  H = H.mpaxis(invpax);
267  H += G;
268  G *= -(2.0*ld+1.0)/ld;
269  G += H;
270 
271  for (size_t i(0); i < G.numx(); ++i) G(0,i) = 0.0;
272  for (size_t i(0); i < H.numx(); ++i) H(0,i) = f(1,i) * Hp0[el];
273 }
274 //--------------------------------------------------------------
275 
276 //--------------------------------------------------------------
277 // Calculation of G00 = -2*Dp* df/dp(p0)
278 void Electric_Field_1D_f1::MakeG00(SHarmonic1D& f) {
279 //--------------------------------------------------------------
280  G = f; G = G.Dp();
281 
282  complex<double> p0p1_sq( pr[0]*pr[0]/(pr[1]*pr[1]) ),
283  inv_mp0p1_sq( 1.0/(1.0-p0p1_sq) ),
284  g_r = -4.0*(pr[1]-pr[0]) * pr[0]/(pr[1]*pr[1]),
285  f00;
286 
287  for (size_t i(0); i < f.numx(); ++i) {
288  f00 = ( f(0,i) - f(1,i) * p0p1_sq) * inv_mp0p1_sq;
289  G(0,i) = ( f(1,i) - f00) * g_r;
290  }
291 }
292 //--------------------------------------------------------------
293 //**************************************************************
294 
295 
296 //**************************************************************
297 //--------------------------------------------------------------
298 Magnetic_Field_1D_f1::Magnetic_Field_1D_f1(size_t Nl, size_t Nm,
299  double pmin, double pmax, size_t Np,
300  double xmin, double xmax, size_t Nx)
301 //--------------------------------------------------------------
302 // Constructor
303 //--------------------------------------------------------------
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 }
343 //--------------------------------------------------------------
344 
345 
346 //--------------------------------------------------------------
347 void Magnetic_Field_1D_f1::operator()(const DistFunc1D& Din,
348  const Field1D& FBx, const Field1D& FBy, const Field1D& FBz,
349  DistFunc1D& Dh) {
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 }
420 //--------------------------------------------------------------
421 //--------------------------------------------------------------
422 void Magnetic_Field_1D_f1::implicit(DistFunc1D& Din,
423  const Field1D& FBx, const Field1D& FBy, const Field1D& FBz,
424  double dt) {
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 }
568 //**************************************************************
569 
570 
571 //**************************************************************
572 //--------------------------------------------------------------
573 Spatial_Advection_1D_f1::Spatial_Advection_1D_f1(size_t Nl, size_t Nm,
574  double pmin, double pmax, size_t Np,
575  double xmin, double xmax, size_t Nx)
576 //--------------------------------------------------------------
577 // Constructor
578 //--------------------------------------------------------------
579  : fd1(Np,Nx),
580  vr(Algorithms::MakeAxis(static_cast<complex<double> >(pmin),
581  static_cast<complex<double> >(pmax),
582  Np)) {
583 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
584 
585  complex<double> lc, mc;
586 
587  for (size_t i(0); i < vr.size(); ++i) {
588  vr[i] = vr[i]/(sqrt(1.0+vr[i]*vr[i]));
589  }
590 
591  double idx = (-1.0) / (2.0*(xmax-xmin)/double(Nx)); // -1/(2dx)
592 
593 // Calculate the "A1, A2" parameters
594  A00 = static_cast<complex<double>>(idx);
595  A10 = static_cast<complex<double>>(idx/3.0);
596  A20 = static_cast<complex<double>>(idx*2.0/5.0);
597 
598 }
599 //--------------------------------------------------------------
600 
601 //--------------------------------------------------------------
602 // Advection in x
603 void Spatial_Advection_1D_f1::operator()(const DistFunc1D& Din, DistFunc1D& Dh) {
604 //--------------------------------------------------------------
605 
606  valarray<complex<double> > vt(vr); vt *= 1.0/Din.mass();
607 
608  fd1 = Din(0,0);
609  fd1 = fd1.Dx(); vt *= A00;
610  Dh(1,0)+=(fd1.mpaxis(vt));
611 
612 
613  fd1 = Din(1,0);
614  fd1 = fd1.Dx(); vt *= A10/A00;
615  Dh(0,0)+=(fd1.mpaxis(vt));
616 
617  fd1 = Din(0,0); fd1 *= static_cast<complex<double>>(1.0/3.0);
618  fd1 = fd1.Dx(); vt *= A20/A10;
619  Dh(1,0)+=(fd1.mpaxis(vt));
620 
621 
622 }
623 //--------------------------------------------------------------
624 
625 
626 //**************************************************************
627 //--------------------------------------------------------------
628 // Functor to be used in the Runge-Kutta methods with explicit
629 // E-field solver
630 
631 //--------------------------------------------------------------
632 // Constructor
633 VlasovFunctor1D_f1_explicitE::VlasovFunctor1D_f1_explicitE(vector<size_t> Nl,vector<size_t> Nm,vector<double> pmax, vector<size_t> Np,
634  double xmin, double xmax, size_t Nx) {
635 //--------------------------------------------------------------
636 
637  for (size_t s(0); s < Nl.size(); ++s){
638 
639  double pmin( pmax[s] / ( double(Np[s] * 2 - 1)) );
640 
641  SA.push_back( Spatial_Advection_1D_f1(Nl[s], Nm[s], pmin, pmax[s], Np[s], xmin, xmax, Nx) );
642 
643  EF.push_back( Electric_Field_1D_f1(Nl[s], Nm[s], pmin, pmax[s], Np[s], xmin, xmax, Nx) );
644 
645  JX.push_back( Current_1D(pmin, pmax[s], Np[s], Nx) );
646 
647  BF.push_back( Magnetic_Field_1D_f1(Nl[s], Nm[s], pmin, pmax[s], Np[s], xmin, xmax, Nx) );
648 
649 // HA.push_back( Hydro_Advection_1D(Nl[s], Nm[s], pmin, pmax[s], Np[s], xmin, xmax, Nx) );
650 
651  AM.push_back( Ampere_1D(xmin, xmax, Nx) );
652 
653  FA.push_back( Faraday_1D(xmin, xmax, Nx) );
654 
655  }
656 }
657 //--------------------------------------------------------------
658 
659 
660 //--------------------------------------------------------------
661 // Collect all of the terms
662 void VlasovFunctor1D_f1_explicitE::operator()(const State1D& Yin, State1D& Yslope){
663 //--------------------------------------------------------------
664 
665  Yslope = 0.0;
666 
667  for (size_t s(0); s < Yin.Species(); ++s) {
668 
669  SA[s](Yin.DF(s),Yslope.DF(s));
670 
671  EF[s](Yin.DF(s),Yin.EMF().Ex(),Yin.EMF().Ey(),Yin.EMF().Ez(),Yslope.DF(s));
672 
673  JX[s](Yin.DF(s),Yslope.EMF().Ex(),Yslope.EMF().Ey(),Yslope.EMF().Ez());
674 
675 // if (Input::List().hydromotion) HA[s](Yin.DF(s),Yin.HYDRO(),Yslope.DF(s));
676 
677  BF[s](Yin.DF(s),Yin.EMF().Bx(),Yin.EMF().By(),Yin.EMF().Bz(),Yslope.DF(s));
678 
679  AM[s](Yin.EMF(),Yslope.EMF());
680 
681  FA[s](Yin.EMF(),Yslope.EMF());
682 
683  Yslope.DF(s) = Yslope.DF(s).Filterp();
684  }
685 
686 }
687 
688 void VlasovFunctor1D_f1_explicitE::operator()(const State1D& Yin, State1D& Yslope, size_t direction){}
689 void VlasovFunctor1D_f1_explicitE::operator()(const State1D& Yin, const State1D& Y2in, State1D& Yslope){}
690 
691 //--------------------------------------------------------------
692 // Constructor
693 VlasovFunctor1D_f1_explicitEB::VlasovFunctor1D_f1_explicitEB(vector<size_t> Nl,vector<size_t> Nm,vector<double> pmax, vector<size_t> Np,
694  double xmin, double xmax, size_t Nx) {
695 //--------------------------------------------------------------
696 
697  for (size_t s(0); s < Nl.size(); ++s){
698 
699  double pmin( pmax[s] / ( double(Np[s] * 2 - 1)) );
700 
701  SA.push_back( Spatial_Advection_1D_f1(Nl[s], Nm[s], pmin, pmax[s], Np[s], xmin, xmax, Nx) );
702 
703  EF.push_back( Electric_Field_1D_f1(Nl[s], Nm[s], pmin, pmax[s], Np[s], xmin, xmax, Nx) );
704 
705  JX.push_back( Current_1D(pmin, pmax[s], Np[s], Nx) );
706 
707 // BF.push_back( Magnetic_Field_1D(Nl[s], Nm[s], pmin, pmax[s], Np[s], xmin, xmax, Nx) );
708 
709 // HA.push_back( Hydro_Advection_1D(Nl[s], Nm[s], pmin, pmax[s], Np[s], xmin, xmax, Nx) );
710 
711  AM.push_back( Ampere_1D(xmin, xmax, Nx) );
712 
713  FA.push_back( Faraday_1D(xmin, xmax, Nx) );
714 
715  }
716 }
717 //--------------------------------------------------------------
718 
719 
720 //--------------------------------------------------------------
721 // Collect all of the terms
722 void VlasovFunctor1D_f1_explicitEB::operator()(const State1D& Yin, State1D& Yslope){
723 //--------------------------------------------------------------
724 
725  Yslope = 0.0;
726 
727  for (size_t s(0); s < Yin.Species(); ++s) {
728 
729  SA[s](Yin.DF(s),Yslope.DF(s));
730 
731  EF[s](Yin.DF(s),Yin.EMF().Ex(),Yin.EMF().Ey(),Yin.EMF().Ez(),Yslope.DF(s));
732 
733  JX[s](Yin.DF(s),Yslope.EMF().Ex(),Yslope.EMF().Ey(),Yslope.EMF().Ez());
734 
735 // if (Input::List().hydromotion) HA[s](Yin.DF(s),Yin.HYDRO(),Yslope.DF(s));
736 
737 // BF[s](Yin.DF(s),Yin.EMF().Bx(),Yin.EMF().By(),Yin.EMF().Bz(),Yslope.DF(s));
738 
739  AM[s](Yin.EMF(),Yslope.EMF());
740 
741  FA[s](Yin.EMF(),Yslope.EMF());
742 
743  // std::cout << "\n";
744 
745  // for (int p(0);p < Yslope.SH(0,1,0).nump(); ++p)
746  // {
747  // std::cout << Yslope.SH(0,1,0)(p,43).real() << "\n";
748  // }
749 
750  // exit(1);
751  Yslope.DF(s) = Yslope.DF(s).Filterp();
752  }
753 
754 }
755 
756 void VlasovFunctor1D_f1_explicitEB::operator()(const State1D& Yin, State1D& Yslope, size_t direction){}
757 void VlasovFunctor1D_f1_explicitEB::operator()(const State1D& Yin, const State1D& Y2in, State1D& Yslope){}
758 //--------------------------------------------------------------
759 // Functor to be used in the Runge-Kutta methods with implicit
760 // E-field solver
761 
762 //--------------------------------------------------------------
763 // Constructor
764 VlasovFunctor1D_f1_implicitE_p1::VlasovFunctor1D_f1_implicitE_p1(vector<size_t> Nl,vector<size_t> Nm,vector<double> pmax, vector<size_t> Np,
765  double xmin, double xmax, size_t Nx) {
766 // //--------------------------------------------------------------
767 
768  for (size_t s(0); s < Nl.size(); ++s){
769 
770  double pmin( pmax[s] / ( double(Np[s] * 2 - 1)) );
771 
772  SA.push_back( Spatial_Advection_1D_f1(Nl[s], Nm[s], pmin, pmax[s], Np[s], xmin, xmax, Nx) );
773 
774  // EF.push_back( Electric_Field_1D(Nl[s], Nm[s], pmin, pmax[s], Np[s], xmin, xmax, Nx) );
775 
776  // JX.push_back( Current_1D(pmin, pmax[s], Np[s], Nx) );
777 
778  // AM.push_back( Ampere_1D(xmin, xmax, Nx) );
779 
780 // HA.push_back( Hydro_Advection_1D(Nl[s], Nm[s], pmin, pmax[s], Np[s], xmin, xmax, Nx) );
781 
782  BF.push_back( Magnetic_Field_1D_f1(Nl[s], Nm[s], pmin, pmax[s], Np[s], xmin, xmax, Nx) );
783 
784 // FA.push_back( Faraday_1D(xmin, xmax, Nx) );
785 
786  }
787 
788 }
789 
790 //--------------------------------------------------------------
791 //
792 //
793 void VlasovFunctor1D_f1_implicitE_p1::operator()(const State1D& Yin, State1D& Yslope){
794 //--------------------------------------------------------------
795 
796  Yslope = 0.0;
797 
798  for (size_t s(0); s < Yin.Species(); ++s) {
799  // Yslope.DF(s).checknan();std::cout << "Vlasov 1 \n";
800  SA[s](Yin.DF(s),Yslope.DF(s));
801  // Yslope.DF(s).checknan();std::cout << "Vlasov 2 \n";
802 
803  // Yslope.DF(s).checknan();std::cout << "Vlasov 3 \n";
804 
805  // EF[s](Yin.DF(s),Yin.EMF().Ex(),Yin.EMF().Ey(),Yin.EMF().Ez(),Yslope.DF(s));
806 
807  // JX[s](Yin.DF(s),Yslope.EMF().Ex(),Yslope.EMF().Ey(),Yslope.EMF().Ez());
808 
809  // AM[s](Yin.EMF(),Yslope.EMF());
810 
811 // if (Input::List().hydromotion) HA[s](Yin.DF(s),Yin.HYDRO(),Yslope.DF(s));
812 
813  BF[s](Yin.DF(s),Yin.EMF().Bx(),Yin.EMF().By(),Yin.EMF().Bz(),Yslope.DF(s));
814  // Yslope.DF(s).checknan();std::cout << "Vlasov 4 \n";
815 // FA[s](Yin.EMF(),Yslope.EMF());
816 
817  // Yslope.DF(s).checknan();std::cout << "Vlasov 5 \n";
818  Yslope.DF(s) = Yslope.DF(s).Filterp();
819 
820 
821  }
822 
823 }
824 void VlasovFunctor1D_f1_implicitE_p1::operator()(const State1D& Yin, State1D& Yslope, size_t direction){}
825 void VlasovFunctor1D_f1_implicitE_p1::operator()(const State1D& Yin, const State1D& Y2in, State1D& Yslope){}
826 //--------------------------------------------------------------
827 // Constructor
828 VlasovFunctor1D_f1_implicitEB_p1::VlasovFunctor1D_f1_implicitEB_p1(vector<size_t> Nl,vector<size_t> Nm,vector<double> pmax, vector<size_t> Np,
829  double xmin, double xmax, size_t Nx) {
830 // //--------------------------------------------------------------
831 
832  for (size_t s(0); s < Nl.size(); ++s){
833 
834  double pmin( pmax[s] / ( double(Np[s] * 2 - 1)) );
835 
836  SA.push_back( Spatial_Advection_1D_f1(Nl[s], Nm[s], pmin, pmax[s], Np[s], xmin, xmax, Nx) );
837 
838  // EF.push_back( Electric_Field_1D(Nl[s], Nm[s], pmin, pmax[s], Np[s], xmin, xmax, Nx) );
839 
840  // JX.push_back( Current_1D(pmin, pmax[s], Np[s], Nx) );
841 
842  // AM.push_back( Ampere_1D(xmin, xmax, Nx) );
843 
844 // HA.push_back( Hydro_Advection_1D(Nl[s], Nm[s], pmin, pmax[s], Np[s], xmin, xmax, Nx) );
845 
846 // BF.push_back( Magnetic_Field_1D(Nl[s], Nm[s], pmin, pmax[s], Np[s], xmin, xmax, Nx) );
847 
848 // FA.push_back( Faraday_1D(xmin, xmax, Nx) );
849 
850  }
851 
852 }
853 
854 //--------------------------------------------------------------
855 //
856 //
857 void VlasovFunctor1D_f1_implicitEB_p1::operator()(const State1D& Yin, State1D& Yslope){
858 //--------------------------------------------------------------
859 
860  Yslope = 0.0;
861 
862  for (size_t s(0); s < Yin.Species(); ++s) {
863  // Yslope.DF(s).checknan();std::cout << "Vlasov 1 \n";
864  SA[s](Yin.DF(s),Yslope.DF(s));
865  // Yslope.DF(s).checknan();std::cout << "Vlasov 2 \n";
866 
867  // Yslope.DF(s).checknan();std::cout << "Vlasov 3 \n";
868 
869  // EF[s](Yin.DF(s),Yin.EMF().Ex(),Yin.EMF().Ey(),Yin.EMF().Ez(),Yslope.DF(s));
870 
871  // JX[s](Yin.DF(s),Yslope.EMF().Ex(),Yslope.EMF().Ey(),Yslope.EMF().Ez());
872 
873  // AM[s](Yin.EMF(),Yslope.EMF());
874 
875 // if (Input::List().hydromotion) HA[s](Yin.DF(s),Yin.HYDRO(),Yslope.DF(s));
876 
877 // BF[s](Yin.DF(s),Yin.EMF().Bx(),Yin.EMF().By(),Yin.EMF().Bz(),Yslope.DF(s));
878  // Yslope.DF(s).checknan();std::cout << "Vlasov 4 \n";
879 // FA[s](Yin.EMF(),Yslope.EMF());
880 
881  // Yslope.DF(s).checknan();std::cout << "Vlasov 5 \n";
882  Yslope.DF(s) = Yslope.DF(s).Filterp();
883 
884 
885  }
886 
887 }
888 void VlasovFunctor1D_f1_implicitEB_p1::operator()(const State1D& Yin, State1D& Yslope, size_t direction){}
889 void VlasovFunctor1D_f1_implicitEB_p1::operator()(const State1D& Yin, const State1D& Y2in, State1D& Yslope){}
890 //--------------------------------------------------------------
891 // Constructor
892 VlasovFunctor1D_f1_implicitE_p2::VlasovFunctor1D_f1_implicitE_p2(vector<size_t> Nl,vector<size_t> Nm,vector<double> pmax, vector<size_t> Np,
893  double xmin, double xmax, size_t Nx) {
894 // //--------------------------------------------------------------
895 
896  for (size_t s(0); s < Nl.size(); ++s){
897 
898  double pmin( pmax[s] / ( double(Np[s] * 2 - 1)) );
899  FA.push_back( Faraday_1D(xmin, xmax, Nx) );
900  EF.push_back( Electric_Field_1D_f1(Nl[s], Nm[s], pmin, pmax[s], Np[s], xmin, xmax, Nx) );
901 
902  }
903 
904 }
905 
906 
907 //------------------------------------------------------------------------------------------------------
908 void VlasovFunctor1D_f1_implicitE_p2::operator()(const State1D& Yin, State1D& Yslope){
909 // //---------------------------------------------------------------------------------------------------
910 
911  Yslope = 0.0;
912 
913  for (size_t s(0); s < Yin.Species(); ++s) {
914  FA[s](Yin.EMF(),Yslope.EMF());
915  EF[s](Yin.DF(s),Yin.EMF().Ex(),Yin.EMF().Ey(),Yin.EMF().Ez(),Yslope.DF(s));
916  Yslope.DF(s) = Yslope.DF(s).Filterp();
917 
918  }
919 
920 }
921 
922 // //--------------------------------------------------------------------------------------------------
923 void VlasovFunctor1D_f1_implicitE_p2::operator()(const State1D& Yin, State1D& Yslope, size_t direction){
924 // //--------------------------------------------------------------------------------------------------
925 
926  Yslope = 0.0;
927 
928  if (direction == 1)
929  {
930  for (size_t s(0); s < Yin.Species(); ++s) {
931  EF[s].Implicit_Ex(Yin.DF(s),Yin.EMF().Ex(),Yslope.DF(s));
932  Yslope.DF(s) = Yslope.DF(s).Filterp();
933  }
934  }
935  else if (direction == 2)
936  {
937  for (size_t s(0); s < Yin.Species(); ++s) {
938  EF[s].Implicit_Ey(Yin.DF(s),Yin.EMF().Ey(),Yslope.DF(s));
939  Yslope.DF(s) = Yslope.DF(s).Filterp();
940  }
941  }
942  else
943  {
944  // complex<double> ii(0.0,1.0);
945  for (size_t s(0); s < Yin.Species(); ++s) {
946  EF[s].Implicit_Ez(Yin.DF(s),Yin.EMF().Ez(),Yslope.DF(s));
947  Yslope.DF(s) = Yslope.DF(s).Filterp();
948  }
949  }
950 
951 }
952 void VlasovFunctor1D_f1_implicitE_p2::operator()(const State1D& Yin, const State1D& Y2in, State1D& Yslope){}
953 // //**************************************************************
Electric_Field_1D_f1::Implicit_Ez
void Implicit_Ez(const DistFunc1D &Din, const Field1D &FEz, DistFunc1D &Dh)
Definition: vlasov_f1.cpp:212
lib-algorithms.h
Algorithms
Electric_Field_1D_f1::pr
valarray< complex< double > > pr
Definition: vlasov_f1.h:69
formulary.h
Plasma Formulary and Units - Declarations.
Magnetic_Field_1D_f1::B1
valarray< complex< double > > B1
Definition: vlasov_f1.h:94
Electric_Field_1D_f1::invpr
valarray< complex< double > > invpr
Definition: vlasov_f1.h:69
Magnetic_Field_1D_f1::implicit
void implicit(DistFunc1D &Din, const Field1D &FBx, const Field1D &FBy, const Field1D &FBz, double dt)
Definition: vlasov_f1.cpp:422
Thomas_Tridiagonal
bool Thomas_Tridiagonal(Array2D< double > &A, valarray< double > &d, valarray< double > &xk)
The tridiagonal solver for implicit collisions.
Definition: nmethods.cpp:216
Magnetic_Field_1D_f1::FLM
SHarmonic1D FLM
Definition: vlasov_f1.h:92
Electric_Field_1D_f1::Electric_Field_1D_f1
Electric_Field_1D_f1(size_t Nl, size_t Nm, double pmin, double pmax, size_t Np, double xmin, double xmax, size_t Nx)
Definition: vlasov_f1.cpp:37
EMF1D::By
Field1D & By()
Definition: state.h:294
Electric_Field_1D_f1::MakeG00
void MakeG00(SHarmonic1D &f)
Definition: vlasov_f1.cpp:278
input.h
Input reader - Declarations.
VlasovFunctor1D_f1_implicitE_p1::VlasovFunctor1D_f1_implicitE_p1
VlasovFunctor1D_f1_implicitE_p1(vector< size_t > Nl, vector< size_t > Nm, vector< double > pmax, vector< size_t > Np, double xmin, double xmax, size_t Nx)
Definition: vlasov_f1.cpp:764
Field1D
A 1D Field.
Definition: state.h:184
lib-array.h
Underlying data structures.
Magnetic_Field_1D_f1::A1
valarray< complex< double > > A1
Definition: vlasov_f1.h:94
Electric_Field_1D_f1::A100
complex< double > A100
Definition: vlasov_f1.h:63
Electric_Field_1D_f1::B211
complex< double > B211
Definition: vlasov_f1.h:63
Electric_Field_1D_f1::Implicit_Ex
void Implicit_Ex(const DistFunc1D &Din, const Field1D &FEx, DistFunc1D &Dh)
Definition: vlasov_f1.cpp:142
DistFunc1D::Filterp
DistFunc1D & Filterp()
Definition: state.cpp:934
Spatial_Advection_1D_f1::Spatial_Advection_1D_f1
Spatial_Advection_1D_f1(size_t Nl, size_t Nm, double pmin, double pmax, size_t Np, double xmin, double xmax, size_t Nx)
Definition: vlasov_f1.cpp:573
VlasovFunctor1D_f1_implicitE_p1::operator()
void operator()(const State1D &Yin, State1D &Yslope)
Definition: vlasov_f1.cpp:793
DistFunc1D::mass
double mass() const
Definition: state.h:400
Magnetic_Field_1D_f1::A3
complex< double > A3
Definition: vlasov_f1.h:96
Electric_Field_1D_f1::C311
complex< double > C311
Definition: vlasov_f1.h:63
SHarmonic1D::mxaxis
SHarmonic1D & mxaxis(const valarray< complex< double > > &shmulti)
Definition: state.cpp:126
Electric_Field_1D_f1::TMP
SHarmonic1D TMP
Definition: vlasov_f1.h:61
SHarmonic1D::Re
SHarmonic1D & Re()
Definition: state.cpp:130
vlasov.h
Functors for various time-integration methodds - Declarations.
SHarmonic1D::numx
size_t numx() const
Definition: state.h:72
VlasovFunctor1D_f1_explicitE::operator()
void operator()(const State1D &Yin, State1D &Yslope)
Definition: vlasov_f1.cpp:662
Electric_Field_1D_f1::MakeGH
void MakeGH(SHarmonic1D &f, size_t l)
Definition: vlasov_f1.cpp:257
Spatial_Advection_1D_f1::A10
complex< double > A10
Definition: vlasov_f1.h:34
EMF1D::Ez
Field1D & Ez()
Definition: state.h:292
SHarmonic1D
A 1D Spherical Harmonic.
Definition: state.h:57
VlasovFunctor1D_f1_explicitEB::operator()
void operator()(const State1D &Yin, State1D &Yslope)
Definition: vlasov_f1.cpp:722
VlasovFunctor1D_f1_implicitE_p2::operator()
void operator()(const State1D &Yin, State1D &Yslope)
Definition: vlasov_f1.cpp:908
SHarmonic1D::Dp
SHarmonic1D & Dp()
Definition: state.cpp:139
Faraday_1D
Definition: vlasov.h:136
Field1D::array
valarray< complex< double > > & array() const
Definition: state.h:196
VlasovFunctor1D_f1_implicitEB_p1::operator()
void operator()(const State1D &Yin, State1D &Yslope)
Definition: vlasov_f1.cpp:857
State1D
Definition: state.h:577
EMF1D::Bx
Field1D & Bx()
Definition: state.h:293
Algorithms::MakeAxis
valarray< T > MakeAxis(const T min, const T max, const size_t N)
Definition: lib-algorithms.h:25
DistFunc1D::m0
size_t m0() const
Definition: state.h:397
DistFunc1D::q
double q() const
Definition: state.h:399
EMF1D::Bz
Field1D & Bz()
Definition: state.h:295
State1D::DF
DistFunc1D & DF(size_t s)
Definition: state.h:602
VlasovFunctor1D_f1_explicitE::VlasovFunctor1D_f1_explicitE
VlasovFunctor1D_f1_explicitE(vector< size_t > Nl, vector< size_t > Nm, vector< double > pmax, vector< size_t > Np, double xmin, double xmax, size_t Nx)
Definition: vlasov_f1.cpp:633
Magnetic_Field_1D_f1
Definition: vlasov_f1.h:75
state.h
Fields, Distributions, Harmonics, States - Declarations.
Spatial_Advection_1D_f1::A20
complex< double > A20
Definition: vlasov_f1.h:34
EMF1D::Ex
Field1D & Ex()
Definition: state.h:290
Spatial_Advection_1D_f1::fd1
SHarmonic1D fd1
Definition: vlasov_f1.h:31
Electric_Field_1D_f1::Implicit_Ey
void Implicit_Ey(const DistFunc1D &Din, const Field1D &FEy, DistFunc1D &Dh)
Definition: vlasov_f1.cpp:171
Electric_Field_1D_f1::G
SHarmonic1D G
Definition: vlasov_f1.h:61
Electric_Field_1D_f1::H
SHarmonic1D H
Definition: vlasov_f1.h:61
Electric_Field_1D_f1::operator()
void operator()(const DistFunc1D &Din, const Field1D &FEx, const Field1D &FEy, const Field1D &FEz, DistFunc1D &Dh)
Definition: vlasov_f1.cpp:81
Electric_Field_1D_f1::C100
complex< double > C100
Definition: vlasov_f1.h:63
SHarmonic1D::Dx
SHarmonic1D & Dx()
Definition: state.cpp:185
VlasovFunctor1D_f1_implicitEB_p1::VlasovFunctor1D_f1_implicitEB_p1
VlasovFunctor1D_f1_implicitEB_p1(vector< size_t > Nl, vector< size_t > Nm, vector< double > pmax, vector< size_t > Np, double xmin, double xmax, size_t Nx)
Definition: vlasov_f1.cpp:828
Magnetic_Field_1D_f1::operator()
void operator()(const DistFunc1D &Din, const Field1D &FBx, const Field1D &FBy, const Field1D &FBz, DistFunc1D &Dh)
Definition: vlasov_f1.cpp:347
EMF1D::Ey
Field1D & Ey()
Definition: state.h:291
Electric_Field_1D_f1::A210
complex< double > A210
Definition: vlasov_f1.h:63
DistFunc1D::l0
size_t l0() const
Definition: state.h:396
Spatial_Advection_1D_f1::operator()
void operator()(const DistFunc1D &Din, DistFunc1D &Dh)
Definition: vlasov_f1.cpp:603
VlasovFunctor1D_f1_explicitEB::VlasovFunctor1D_f1_explicitEB
VlasovFunctor1D_f1_explicitEB(vector< size_t > Nl, vector< size_t > Nm, vector< double > pmax, vector< size_t > Np, double xmin, double xmax, size_t Nx)
Definition: vlasov_f1.cpp:693
Spatial_Advection_1D_f1::A00
complex< double > A00
Definition: vlasov_f1.h:34
Electric_Field_1D_f1
Definition: vlasov_f1.h:42
Electric_Field_1D_f1::Hp0
valarray< complex< double > > Hp0
Definition: vlasov_f1.h:69
Spatial_Advection_1D_f1::vr
valarray< complex< double > > vr
Definition: vlasov_f1.h:33
Electric_Field_1D_f1::A310
complex< double > A310
Definition: vlasov_f1.h:63
State1D::Species
size_t Species() const
Definition: state.h:596
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: vlasov_f1.cpp:298
DistFunc1D
The 1D distribution function is the container for all SHarmonic1D per species.
Definition: state.h:376
vlasov_f1.h
Vlasov Equation - Declarations.
Current_1D
Definition: vlasov.h:115
Ampere_1D
Definition: vlasov.h:154
State1D::EMF
EMF1D & EMF() const
Definition: state.h:610
VlasovFunctor1D_f1_implicitE_p2::VlasovFunctor1D_f1_implicitE_p2
VlasovFunctor1D_f1_implicitE_p2(vector< size_t > Nl, vector< size_t > Nm, vector< double > pmax, vector< size_t > Np, double xmin, double xmax, size_t Nx)
Definition: vlasov_f1.cpp:892
SHarmonic1D::mpaxis
SHarmonic1D & mpaxis(const valarray< complex< double > > &shmulti)
Definition: state.cpp:122
nmethods.h
Numerical Methods - Declarations.