OSHUN  beta
Arbitrary Order Spherical-Harmonic 1D-3P Vlasov-Fokker-Planck-Maxwell code
functors.cpp
Go to the documentation of this file.
1 
8 //--------------------------------------------------------------
9 // Standard libraries
10 #include <iostream>
11 #include <vector>
12 #include <valarray>
13 #include <complex>
14 #include <algorithm>
15 #include <cstdlib>
16 
17 #include <math.h>
18 #include <map>
19 
20 // My libraries
21 #include "lib-array.h"
22 #include "lib-algorithms.h"
23 #include "nmethods.h"
24 
25 // Declerations
26 #include "state.h"
27 #include "formulary.h"
28 #include "input.h"
29 #include "fluid.h"
30 #include "vlasov.h"
31 #include "functors.h"
32 
33 //**************************************************************
34 
35 //**************************************************************
36 //--------------------------------------------------------------
37 // Functor to be used in the Runge-Kutta methods with explicit
38 // E-field solver
39 
40 //--------------------------------------------------------------
41 // Constructor
42 VlasovFunctor1D_explicitE::VlasovFunctor1D_explicitE(vector<size_t> Nl,vector<size_t> Nm,vector<double> pmax, vector<size_t> Np,
43  double xmin, double xmax, size_t Nx) {
44 //--------------------------------------------------------------
45 
46  for (size_t s(0); s < Nl.size(); ++s){
47 
48  SA.push_back( Spatial_Advection_1D(Nl[s], Nm[s], 0.0, pmax[s], Np[s], xmin, xmax, Nx) );
49 
50  EF.push_back( Electric_Field_1D(Nl[s], Nm[s], 0.0, pmax[s], Np[s], xmin, xmax, Nx) );
51 
52  JX.push_back( Current_1D(0.0, pmax[s], Np[s], Nx) );
53 
54  BF.push_back( Magnetic_Field_1D(Nl[s], Nm[s], 0.0, pmax[s], Np[s], xmin, xmax, Nx) );
55 
56  AM.push_back( Ampere_1D(xmin, xmax, Nx) );
57 
58  FA.push_back( Faraday_1D(xmin, xmax, Nx) );
59 
60  }
61 }
62 //--------------------------------------------------------------
63 
64 
65 //--------------------------------------------------------------
66 // Collect all of the terms
68 //--------------------------------------------------------------
69  bool debug(0);
70 
71  Yslope = 0.0;
72 
73  for (size_t s(0); s < Yin.Species(); ++s) {
74 
75  if (Yin.DF(s).l0() == 1) {
76 
77  SA[s].f1only(Yin.DF(s),Yslope.DF(s));
78 
79  EF[s].f1only(Yin.DF(s),Yin.EMF().Ex(),Yin.EMF().Ey(),Yin.EMF().Ez(),Yslope.DF(s));
80 
81  BF[s].f1only(Yin.DF(s),Yin.EMF().Bx(),Yin.EMF().By(),Yin.EMF().Bz(),Yslope.DF(s));
82 
83  JX[s](Yin.DF(s),Yslope.EMF().Ex(),Yslope.EMF().Ey(),Yslope.EMF().Ez());
84 
85  AM[s](Yin.EMF(),Yslope.EMF());
86 
87  FA[s](Yin.EMF(),Yslope.EMF());
88 
89  }
90  else if (Yin.DF(s).m0() == 0) {
91 
92 
93  if (debug)
94  {
95  std::cout << "\n\n f at start:";
96  for (size_t ip(0); ip < Yin.SH(0,0,0).nump(); ++ip){
97  std::cout << "\nf(" << ip << ") = " << Yin.SH(0,1,0)(ip,4);
98  }
99 
100  std::cout << "\n\n E at start:";
101  for (size_t ix(0); ix < Yin.SH(0,0,0).numx(); ++ix){
102  std::cout << "\nEx(" << ix << ") = " << Yin.EMF().Ex()(ix);
103  }
104  }
105  EF[s].es1d(Yin.DF(s),Yin.EMF().Ex(),Yin.EMF().Ey(),Yin.EMF().Ez(),Yslope.DF(s));
106 
107 
108  if (debug)
109  {
110  std::cout << "\n\nf after E:";
111  for (size_t ip(0); ip < Yin.SH(0,0,0).nump(); ++ip){
112  std::cout << "\nf(" << ip << ") = " << Yslope.SH(0,1,0)(ip,4);
113  }
114  }
115  JX[s].es1d(Yin.DF(s),Yslope.EMF().Ex());
116 
117  if (debug)
118  {
119  std::cout << "\n\n after J:";
120  for (size_t ix(0); ix < Yin.SH(0,0,0).numx(); ++ix){
121  std::cout << "\nEx(" << ix << ") = " << Yslope.EMF().Ex()(ix);
122  }
123  }
124 
125  SA[s].es1d(Yin.DF(s),Yslope.DF(s));
126 
127  if (debug)
128  {
129  std::cout << "\n\n after SA:";
130  for (size_t ip(0); ip < Yin.SH(0,0,0).nump(); ++ip){
131  std::cout << "\nf(" << ip << ") = " << Yslope.SH(0,1,0)(ip,4);
132  }
133  }
134  }
135 
136  else {
137 
138  SA[s](Yin.DF(s),Yslope.DF(s));
139 
140  EF[s](Yin.DF(s),Yin.EMF().Ex(),Yin.EMF().Ey(),Yin.EMF().Ez(),Yslope.DF(s));
141 
142  BF[s](Yin.DF(s),Yin.EMF().Bx(),Yin.EMF().By(),Yin.EMF().Bz(),Yslope.DF(s));
143 
144  JX[s](Yin.DF(s),Yslope.EMF().Ex(),Yslope.EMF().Ey(),Yslope.EMF().Ez());
145 
146  AM[s](Yin.EMF(),Yslope.EMF());
147 
148  FA[s](Yin.EMF(),Yslope.EMF());
149 
150  }
151 
152  // Yslope.DF(s) = Yslope.DF(s).Filterp();
153 
154  // exit(1);
155  }
156 
157 }
158 
159 void VlasovFunctor1D_explicitE::operator()(const State1D& Yin, State1D& Yslope, size_t direction){}
160 //--------------------------------------------------------------
161 // Functors to be used in the Leapfrog methods with explicit
162 // E-field solver
163 
164 //--------------------------------------------------------------
165 // Constructor
166 VlasovFunctor1D_spatialpush::VlasovFunctor1D_spatialpush(vector<size_t> Nl,vector<size_t> Nm,vector<double> pmax, vector<size_t> Np,
167  double xmin, double xmax, size_t Nx) {
168 //--------------------------------------------------------------
169 
170  for (size_t s(0); s < Nl.size(); ++s){
171 
172  SA.push_back( Spatial_Advection_1D(Nl[s], Nm[s], 0.0, pmax[s], Np[s], xmin, xmax, Nx) );
173 
174  }
175 }
176 //--------------------------------------------------------------
177 
178 
179 //--------------------------------------------------------------
180 // Collect all of the terms
182 //--------------------------------------------------------------
183 
184  Yslope = 0.0;
185 
186  for (size_t s(0); s < Yin.Species(); ++s) {
187 
188  if (Yin.DF(s).l0() == 1) {
189 
190  SA[s].f1only(Yin.DF(s),Yslope.DF(s));
191 
192  }
193  else if (Yin.DF(s).m0() == 0) {
194 
195  SA[s].es1d(Yin.DF(s),Yslope.DF(s));
196 
197  }
198 
199  else {
200 
201  SA[s](Yin.DF(s),Yslope.DF(s));
202 
203  }
204 
205  Yslope.DF(s) = Yslope.DF(s).Filterp();
206  }
207 
208 }
209 
210 void VlasovFunctor1D_spatialpush::operator()(const State1D& Yin, State1D& Yslope, size_t direction){}
211 //--------------------------------------------------------------
212 // Constructor
213 VlasovFunctor1D_momentumpush::VlasovFunctor1D_momentumpush(vector<size_t> Nl,vector<size_t> Nm,vector<double> pmax, vector<size_t> Np,
214  double xmin, double xmax, size_t Nx) {
215 //--------------------------------------------------------------
216 
217  for (size_t s(0); s < Nl.size(); ++s){
218 
219  EF.push_back( Electric_Field_1D(Nl[s], Nm[s], 0.0, pmax[s], Np[s], xmin, xmax, Nx) );
220 
221  JX.push_back( Current_1D(0.0, pmax[s], Np[s], Nx) );
222 
223  BF.push_back( Magnetic_Field_1D(Nl[s], Nm[s], 0.0, pmax[s], Np[s], xmin, xmax, Nx) );
224 
225  AM.push_back( Ampere_1D(xmin, xmax, Nx) );
226 
227  FA.push_back( Faraday_1D(xmin, xmax, Nx) );
228 
229  }
230 }
231 //--------------------------------------------------------------
232 
233 
234 //--------------------------------------------------------------
235 // Collect all of the terms
237 //--------------------------------------------------------------
238 
239  Yslope = 0.0;
240 
241  for (size_t s(0); s < Yin.Species(); ++s) {
242 
243  if (Yin.DF(s).l0() == 1) {
244 
245  EF[s].f1only(Yin.DF(s),Yin.EMF().Ex(),Yin.EMF().Ey(),Yin.EMF().Ez(),Yslope.DF(s));
246 
247  BF[s].f1only(Yin.DF(s),Yin.EMF().Bx(),Yin.EMF().By(),Yin.EMF().Bz(),Yslope.DF(s));
248 
249  JX[s](Yin.DF(s),Yslope.EMF().Ex(),Yslope.EMF().Ey(),Yslope.EMF().Ez());
250 
251  AM[s](Yin.EMF(),Yslope.EMF());
252 
253  FA[s](Yin.EMF(),Yslope.EMF());
254 
255 
256 
257  }
258  else if (Yin.DF(s).m0() == 0) {
259 
260  EF[s].es1d(Yin.DF(s),Yin.EMF().Ex(),Yin.EMF().Ey(),Yin.EMF().Ez(),Yslope.DF(s));
261 
262  JX[s].es1d(Yin.DF(s),Yslope.EMF().Ex());
263 
264  }
265 
266  else {
267 
268  EF[s](Yin.DF(s),Yin.EMF().Ex(),Yin.EMF().Ey(),Yin.EMF().Ez(),Yslope.DF(s));
269 
270  BF[s](Yin.DF(s),Yin.EMF().Bx(),Yin.EMF().By(),Yin.EMF().Bz(),Yslope.DF(s));
271 
272  JX[s](Yin.DF(s),Yslope.EMF().Ex(),Yslope.EMF().Ey(),Yslope.EMF().Ez());
273 
274  AM[s](Yin.EMF(),Yslope.EMF());
275 
276  FA[s](Yin.EMF(),Yslope.EMF());
277 
278  }
279 
280  Yslope.DF(s) = Yslope.DF(s).Filterp();
281  }
282 
283 }
284 
285 void VlasovFunctor1D_momentumpush::operator()(const State1D& Yin, State1D& Yslope, size_t direction){}
286 //--------------------------------------------------------------
287 // Constructor
288 VlasovFunctor1D_explicitE_implicitB::VlasovFunctor1D_explicitE_implicitB(vector<size_t> Nl,vector<size_t> Nm,vector<double> pmax, vector<size_t> Np,
289  double xmin, double xmax, size_t Nx) {
290 //--------------------------------------------------------------
291 
292  for (size_t s(0); s < Nl.size(); ++s){
293 
294  SA.push_back( Spatial_Advection_1D(Nl[s], Nm[s], 0.0, pmax[s], Np[s], xmin, xmax, Nx) );
295 
296  EF.push_back( Electric_Field_1D(Nl[s], Nm[s], 0.0, pmax[s], Np[s], xmin, xmax, Nx) );
297 
298  JX.push_back( Current_1D(0.0, pmax[s], Np[s], Nx) );
299 
300  AM.push_back( Ampere_1D(xmin, xmax, Nx) );
301 
302  FA.push_back( Faraday_1D(xmin, xmax, Nx) );
303 
304  }
305 }
306 //--------------------------------------------------------------
307 
308 
309 //--------------------------------------------------------------
310 // Collect all of the terms
312 //--------------------------------------------------------------
313 
314  Yslope = 0.0;
315 
316  for (size_t s(0); s < Yin.Species(); ++s) {
317 
318  if (Yin.DF(s).l0() == 1) {
319  SA[s].f1only(Yin.DF(s),Yslope.DF(s));
320  EF[s].f1only(Yin.DF(s),Yin.EMF().Ex(),Yin.EMF().Ey(),Yin.EMF().Ez(),Yslope.DF(s));
321  }
322  else{
323  SA[s](Yin.DF(s),Yslope.DF(s));
324  EF[s](Yin.DF(s),Yin.EMF().Ex(),Yin.EMF().Ey(),Yin.EMF().Ez(),Yslope.DF(s));
325  }
326 
327  JX[s](Yin.DF(s),Yslope.EMF().Ex(),Yslope.EMF().Ey(),Yslope.EMF().Ez());
328 
329  AM[s](Yin.EMF(),Yslope.EMF());
330 
331  FA[s](Yin.EMF(),Yslope.EMF());
332 
333  Yslope.DF(s) = Yslope.DF(s).Filterp();
334  }
335 
336 }
337 
338 void VlasovFunctor1D_explicitE_implicitB::operator()(const State1D& Yin, State1D& Yslope, size_t direction){}
339 //--------------------------------------------------------------
340 // Functor to be used in the Runge-Kutta methods with implicit
341 // E-field solver
342 
343 //--------------------------------------------------------------
344 // Constructor
345 VlasovFunctor1D_implicitE_p1::VlasovFunctor1D_implicitE_p1(vector<size_t> Nl,vector<size_t> Nm,vector<double> pmax, vector<size_t> Np,
346  double xmin, double xmax, size_t Nx) {
347 // //--------------------------------------------------------------
348 
349  for (size_t s(0); s < Nl.size(); ++s){
350 
351  SA.push_back( Spatial_Advection_1D(Nl[s], Nm[s], 0.0, pmax[s], Np[s], xmin, xmax, Nx) );
352 
353  HA.push_back( Hydro_Advection_1D(Nl[s], Nm[s], 0.0, pmax[s], Np[s], xmin, xmax, Nx) );
354 
355  BF.push_back( Magnetic_Field_1D(Nl[s], Nm[s], 0.0, pmax[s], Np[s], xmin, xmax, Nx) );
356 
357  }
358 
359 }
360 
361 //--------------------------------------------------------------
362 //
363 //
365 //--------------------------------------------------------------
366 
367  Yslope = 0.0;
368 
369  for (size_t s(0); s < Yin.Species(); ++s) {
370 
371 
372  if (Yin.DF(s).l0() == 1) {
373  SA[s].f1only(Yin.DF(s),Yslope.DF(s));
374  BF[s].f1only(Yin.DF(s),Yin.EMF().Bx(),Yin.EMF().By(),Yin.EMF().Bz(),Yslope.DF(s));
375  // HA[s](Yin.DF(s),Yin.HYDRO(),Yslope.DF(s));
376  }
377  else {
378  SA[s](Yin.DF(s),Yslope.DF(s));
379  BF[s](Yin.DF(s),Yin.EMF().Bx(),Yin.EMF().By(),Yin.EMF().Bz(),Yslope.DF(s));
380  // HA[s](Yin.DF(s),Yin.HYDRO(),Yslope.DF(s));
381  }
382 
383  Yslope.DF(s) = Yslope.DF(s).Filterp();
384 
385 
386  }
387 
388 }
389 void VlasovFunctor1D_implicitE_p1::operator()(const State1D& Yin, State1D& Yslope, size_t direction){}
390 //--------------------------------------------------------------
391 // Constructor
392 VlasovFunctor1D_implicitE_implicitB_p1::VlasovFunctor1D_implicitE_implicitB_p1(vector<size_t> Nl,vector<size_t> Nm,vector<double> pmax, vector<size_t> Np,
393  double xmin, double xmax, size_t Nx) {
394 // //--------------------------------------------------------------
395 
396  for (size_t s(0); s < Nl.size(); ++s){
397 
398  SA.push_back( Spatial_Advection_1D(Nl[s], Nm[s], 0.0, pmax[s], Np[s], xmin, xmax, Nx) );
399 
400  HA.push_back( Hydro_Advection_1D(Nl[s], Nm[s], 0.0, pmax[s], Np[s], xmin, xmax, Nx) );
401  }
402 
403 }
404 
405 //--------------------------------------------------------------
406 //
407 //
409 //--------------------------------------------------------------
410 
411  Yslope = 0.0;
412 
413  for (size_t s(0); s < Yin.Species(); ++s) {
414 
415  if (Yin.DF(s).l0() == 1) SA[s].f1only(Yin.DF(s),Yslope.DF(s));
416  else
417  {
418  SA[s](Yin.DF(s),Yslope.DF(s));
419  // HA[s](Yin.DF(s),Yin.HYDRO(),Yslope.DF(s));
420  }
421 
422  Yslope.DF(s) = Yslope.DF(s).Filterp();
423 
424  }
425 
426 }
427 void VlasovFunctor1D_implicitE_implicitB_p1::operator()(const State1D& Yin, State1D& Yslope, size_t direction){}
428 //--------------------------------------------------------------
429 // Constructor
430 VlasovFunctor1D_implicitE_p2::VlasovFunctor1D_implicitE_p2(vector<size_t> Nl,vector<size_t> Nm,vector<double> pmax, vector<size_t> Np,
431  double xmin, double xmax, size_t Nx) {
432 // //--------------------------------------------------------------
433 
434  for (size_t s(0); s < Nl.size(); ++s){
435 
436  FA.push_back( Faraday_1D(xmin, xmax, Nx) );
437  EF.push_back( Electric_Field_1D(Nl[s], Nm[s], 0.0, pmax[s], Np[s], xmin, xmax, Nx) );
438 
439  }
440 
441 }
442 //------------------------------------------------------------------------------------------------------
444 // //---------------------------------------------------------------------------------------------------
445 
446  Yslope = 0.0;
447 
448  for (size_t s(0); s < Yin.Species(); ++s) {
449  FA[s](Yin.EMF(),Yslope.EMF());
450 
451  if (Yin.DF(s).l0() == 1) EF[s].f1only(Yin.DF(s),Yin.EMF().Ex(),Yin.EMF().Ey(),Yin.EMF().Ez(),Yslope.DF(s));
452  else EF[s](Yin.DF(s),Yin.EMF().Ex(),Yin.EMF().Ey(),Yin.EMF().Ez(),Yslope.DF(s));
453 
454  Yslope.DF(s) = Yslope.DF(s).Filterp();
455 
456  }
457 
458 }
459 
460 //--------------------------------------------------------------------------------------------------
461 void VlasovFunctor1D_implicitE_p2::operator()(const State1D& Yin, State1D& Yslope, size_t direction){
462 //--------------------------------------------------------------------------------------------------
463 
464  Yslope = 0.0;
465 
466  if (direction == 1)
467  {
468  for (size_t s(0); s < Yin.Species(); ++s) {
469  if (Yin.DF(s).l0() == 1) EF[s].Implicit_Ex_f1only(Yin.DF(s),Yin.EMF().Ex(),Yslope.DF(s));
470  else EF[s].Implicit_Ex(Yin.DF(s),Yin.EMF().Ex(),Yslope.DF(s));
471  Yslope.DF(s) = Yslope.DF(s).Filterp();
472  }
473  }
474  else if (direction == 2)
475  {
476  for (size_t s(0); s < Yin.Species(); ++s) {
477  if (Yin.DF(s).l0() == 1) EF[s].Implicit_Ey_f1only(Yin.DF(s),Yin.EMF().Ey(),Yslope.DF(s));
478  else EF[s].Implicit_Ey(Yin.DF(s),Yin.EMF().Ey(),Yslope.DF(s));
479  Yslope.DF(s) = Yslope.DF(s).Filterp();
480  }
481  }
482  else
483  {
484  // complex<double> ii(0.0,1.0);
485  for (size_t s(0); s < Yin.Species(); ++s) {
486  if (Yin.DF(s).l0() == 1) EF[s].Implicit_Ez_f1only(Yin.DF(s),Yin.EMF().Ez(),Yslope.DF(s));
487  else EF[s].Implicit_Ez(Yin.DF(s),Yin.EMF().Ez(),Yslope.DF(s));
488  Yslope.DF(s) = Yslope.DF(s).Filterp();
489  }
490  }
491 
492 }
493 //**************************************************************
Algorithms
Plasma Formulary and Units - Declarations.
Field1D & By()
Definition: state.h:294
Input reader - Declarations.
Underlying data structures.
VlasovFunctor1D_spatialpush(vector< size_t > Nl, vector< size_t > Nm, vector< double > pmax, vector< size_t > Np, double xmin, double xmax, size_t Nx)
Definition: functors.cpp:166
void operator()(const State1D &Yin, State1D &Yslope)
Definition: functors.cpp:443
DistFunc1D & Filterp()
Definition: state.cpp:934
VlasovFunctor1D_explicitE_implicitB(vector< size_t > Nl, vector< size_t > Nm, vector< double > pmax, vector< size_t > Np, double xmin, double xmax, size_t Nx)
Definition: functors.cpp:288
VlasovFunctor1D_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: functors.cpp:430
void operator()(const State1D &Yin, State1D &Yslope)
Definition: functors.cpp:67
Functors for various time-integration methodds - Declarations.
size_t numx() const
Definition: state.h:72
Field1D & Ez()
Definition: state.h:292
vector< Faraday_1D > FA
Definition: functors.h:40
Definition: state.h:577
Field1D & Bx()
Definition: state.h:293
VlasovFunctor1D_explicitE(vector< size_t > Nl, vector< size_t > Nm, vector< double > pmax, vector< size_t > Np, double xmin, double xmax, size_t Nx)
Definition: functors.cpp:42
size_t nump() const
Definition: state.h:71
size_t m0() const
Definition: state.h:397
void operator()(const State1D &Yin, State1D &Yslope)
Definition: functors.cpp:408
Field1D & Bz()
Definition: state.h:295
DistFunc1D & DF(size_t s)
Definition: state.h:602
void operator()(const State1D &Yin, State1D &Yslope)
Definition: functors.cpp:311
void operator()(const State1D &Yin, State1D &Yslope)
Definition: functors.cpp:236
SHarmonic1D & SH(size_t s, size_t lh, size_t mh)
Definition: state.h:605
Fields, Distributions, Harmonics, States - Declarations.
Field1D & Ex()
Definition: state.h:290
void operator()(const State1D &Yin, State1D &Yslope)
Definition: functors.cpp:364
VlasovFunctor1D_implicitE_implicitB_p1(vector< size_t > Nl, vector< size_t > Nm, vector< double > pmax, vector< size_t > Np, double xmin, double xmax, size_t Nx)
Definition: functors.cpp:392
vector< Magnetic_Field_1D > BF
Definition: functors.h:39
VlasovFunctor1D_momentumpush(vector< size_t > Nl, vector< size_t > Nm, vector< double > pmax, vector< size_t > Np, double xmin, double xmax, size_t Nx)
Definition: functors.cpp:213
vector< Ampere_1D > AM
Definition: functors.h:38
vector< Current_1D > JX
Definition: functors.h:37
VlasovFunctor1D_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: functors.cpp:345
vector< Spatial_Advection_1D > SA
Definition: functors.h:35
Field1D & Ey()
Definition: state.h:291
size_t l0() const
Definition: state.h:396
void operator()(const State1D &Yin, State1D &Yslope)
Definition: functors.cpp:181
size_t Species() const
Definition: state.h:596
vector< Electric_Field_1D > EF
Definition: functors.h:36
EMF1D & EMF() const
Definition: state.h:610
Numerical Methods - Declarations.