OSHUN  beta
Arbitrary Order Spherical-Harmonic 1D-3P Vlasov-Fokker-Planck-Maxwell code
Functions
nmethods.h File Reference

Numerical Methods - Declarations. More...

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Functions

bool Gauss_Seidel (Array2D< double > &A, valarray< complex< double > > &b, valarray< complex< double > > &xk)
 Performs Gauss-Seidel method on Ax = b. More...
 
void TridiagonalSolve (const valarray< double > &a, const valarray< double > &b, valarray< double > &c, valarray< complex< double > > d, valarray< complex< double > > &x)
 Set up tridiagonal solver. More...
 
void TridiagonalSolve (const valarray< double > &a, const valarray< double > &b, valarray< double > &c, valarray< double > d, valarray< double > &x)
 
void TridiagonalSolve (const valarray< complex< double > > &a, const valarray< complex< double > > &b, valarray< complex< double > > &c, valarray< complex< double > > d, valarray< complex< double > > &x)
 
bool Thomas_Tridiagonal (Array2D< double > &A, valarray< double > &d, valarray< double > &xk)
 The tridiagonal solver for implicit collisions. More...
 
bool Thomas_Tridiagonal (Array2D< double > &A, valarray< complex< double > > &d, valarray< complex< double > > &xk)
 
bool Thomas_Tridiagonal (Array2D< complex< double > > &A, valarray< complex< double > > &d, valarray< complex< double > > &xk)
 
complex< double > Det33 (Array2D< complex< double > > &A)
 
complex< double > Detx33 (valarray< complex< double > > &D, Array2D< complex< double > > &A)
 
complex< double > Dety33 (valarray< complex< double > > &D, Array2D< complex< double > > &A)
 
complex< double > Detz33 (valarray< complex< double > > &D, Array2D< complex< double > > &A)
 
valarray< float > vfloat (const valarray< double > &vDouble)
 
vector< float > vfloat (const vector< double > vDouble)
 Convert data structure to float structure. More...
 
vector< float > vfloat (const vector< complex< double > > vDouble)
 
vector< float > vfloat_complex (const vector< complex< double > > vDouble)
 
valarray< double > df_4thorder (const valarray< double > &f)
 
valarray< double > df_4thorder (valarray< double > &f)
 
Array2D< complex< double > > df1_4thorder (Array2D< complex< double > > &f)
 
Array2D< complex< double > > df2_4thorder (Array2D< complex< double > > &f)
 

Detailed Description

Numerical Methods - Declarations.

Author
PICKSC
Date
September 1, 2016

This cpp file contains the definitions for the functions required for the numerical methods.

Definition in file nmethods.h.

Function Documentation

◆ Det33()

complex<double> Det33 ( Array2D< complex< double > > &  A)

◆ Detx33()

complex<double> Detx33 ( valarray< complex< double > > &  D,
Array2D< complex< double > > &  A 
)

◆ Dety33()

complex<double> Dety33 ( valarray< complex< double > > &  D,
Array2D< complex< double > > &  A 
)

◆ Detz33()

complex<double> Detz33 ( valarray< complex< double > > &  D,
Array2D< complex< double > > &  A 
)

◆ df1_4thorder()

Array2D<complex<double> > df1_4thorder ( Array2D< complex< double > > &  f)

Definition at line 469 of file nmethods.cpp.

469  {
470  Array2D<complex<double> > df(f.dim1(),f.dim2());
471 
472  for (long i2(0); i2<f.dim1();++i2){
473 
474  df(0,i2) = f(1,i2)-f(0,i2);
475  df(1,i2) = 1.0/12.0*(f(4,i2)-6.0*f(3,i2)+18.0*f(2,i2)-10.0*f(1,i2)-3.0*f(0,i2));
476 
477  for (long i1(2); i1<f.dim1()-2;++i1){
478  df(i1,i2) = 1.0/12.0*(-f(i1+2,i2)+8.0*f(i1+1,i2)-8.0*f(i1-1,i2)+f(i1-2,i2));
479  }
480 
481  df(f.dim1()-2,i2) = 1.0/12.0*(3.0*f(f.dim2()-1,i2)+10.0*f(f.dim2()-2,i2)-18.0*f(f.dim2()-3,i2)+6.0*f(f.dim2()-4,i2)-f(f.dim2()-5,i2));
482  df(f.dim1()-1,i2) = f(f.dim2()-1,i2)-f(f.dim2()-2,i2);
483 
484  }
485 
486  return df;
487 }
Array2D::dim1
size_t dim1() const
Definition: lib-array.h:289
Array2D::dim2
size_t dim2() const
Definition: lib-array.h:290

◆ df2_4thorder()

Array2D<complex<double> > df2_4thorder ( Array2D< complex< double > > &  f)

Definition at line 489 of file nmethods.cpp.

489  {
490  Array2D<complex<double> > df(f.dim1(),f.dim2());
491 
492  for (long i1(0); i1<f.dim1();++i1){
493 
494  df(i1,0) = f(i1,1)-f(i1,0);
495  df(i1,1) = 1.0/12.0*(f(i1,4)-6.0*f(i1,3)+18.0*f(i1,2)-10.0*f(i1,1)-3.0*f(i1,0));
496 
497  for (long i2(2); i2<f.dim2()-2;++i2){
498  df(i1,i2) = 1.0/12.0*(-f(i1,i2+2)+8.0*f(i1,i2+1)-8.0*f(i1,i2-1)+f(i1,i2-2));
499  }
500 
501  df(i1,f.dim2()-2) = 1.0/12.0*(3.0*f(i1,f.dim2()-1)+10.0*f(i1,f.dim2()-2)-18.0*f(i1,f.dim2()-3)+6.0*f(i1,f.dim2()-4)-f(i1,f.dim2()-5));
502  df(i1,f.dim2()-1) = f(i1,f.dim2()-1)-f(i1,f.dim2()-2);
503 
504  }
505 
506  return df;
507 }
Array2D::dim1
size_t dim1() const
Definition: lib-array.h:289
Array2D::dim2
size_t dim2() const
Definition: lib-array.h:290

◆ df_4thorder() [1/2]

valarray<double> df_4thorder ( const valarray< double > &  f)

Definition at line 421 of file nmethods.cpp.

Referenced by Fluid_Equation_1D::chargefraction(), Fluid_Equation_1D::density(), and Fluid_Equation_1D::velocity().

421  {
422  valarray<double> df(f.size());
423 
424  df[0] = f[1]-f[0];
425  df[1] = 1.0/12.0*(f[4]-6.0*f[3]+18.0*f[2]-10.0*f[1]-3.0*f[0]);
426 
427  for (size_t i(2); i < df.size()-2; ++i) {
428  df[i] = 1.0/12.0*(-f[i+2]+8.0*f[i+1]-8.0*f[i-1]+f[i-2]);
429  }
430 
431  df[df.size()-2] = 1.0/12.0*(3.0*f[df.size()-1]+10.0*f[df.size()-2]-18.0*f[df.size()-3]+6.0*f[df.size()-4]-f[df.size()-5]);
432  df[df.size()-1] = f[df.size()-1]-f[df.size()-2];
433 
434  return df;
435 }
Here is the caller graph for this function:

◆ df_4thorder() [2/2]

valarray<double> df_4thorder ( valarray< double > &  f)

Definition at line 437 of file nmethods.cpp.

437  {
438  valarray<double> df(f.size());
439 
440  df[0] = f[1]-f[0];
441  df[1] = 1.0/12.0*(f[4]-6.0*f[3]+18.0*f[2]-10.0*f[1]-3.0*f[0]);
442 
443  for (size_t i(2); i < df.size()-2; ++i) {
444  df[i] = 1.0/12.0*(-f[i+2]+8.0*f[i+1]-8.0*f[i-1]+f[i-2]);
445  }
446 
447  df[df.size()-2] = 1.0/12.0*(3.0*f[df.size()-1]+10.0*f[df.size()-2]-18.0*f[df.size()-3]+6.0*f[df.size()-4]-f[df.size()-5]);
448  df[df.size()-1] = f[df.size()-1]-f[df.size()-2];
449 
450  return df;
451 }

◆ Gauss_Seidel()

bool Gauss_Seidel ( Array2D< double > &  A,
valarray< complex< double > > &  b,
valarray< complex< double > > &  xk 
)

Performs Gauss-Seidel method on Ax = b.

Parameters
AInput matrix that contains collisional integrals and coefficients
bVector representing distribution function at current time-step
xkThe solution vector
Returns
The solution vector representing the distribution function at the next time-step

Fills solution into xk. The other matrices are not modified The function returns "false" if the matrix A is not diagonally dominant

Definition at line 29 of file nmethods.cpp.

References Array2D< T >::dim1(), and Array2D< T >::dim2().

Referenced by interspecies_flm_implicit_step::advance(), and self_flm_implicit_step::advance().

31  {
32 //-------------------------------------------------------------------
33  double tol(1.0e-1); //< Tolerance for absolute error
34  int MAXiter(5); //< Maximum iteration allowed
35 
36 
37 // The Matrices all have the right dimensions
38 // -------------------------------------------------------------
39  if ( ( A.dim1() != A.dim2() ) ||
40  ( A.dim1() != b.size() ) ||
41  ( A.dim1() != xk.size() ) ) {
42  cout << "Error: The Matrices don't have the right dimensions!" << endl;
43  exit(1);
44  }
45 // Check if the matrix A is diagonally dominant
46 // -------------------------------------------------------------
47  for (int i(0); i < A.dim1(); ++i){
48  double rowi(0.0);
49  for (int j(0); j < A.dim2(); ++j){
50  rowi += A(i,j);
51  }
52  if (!(rowi < 2.0*A(i,i))) return false;
53  }
54 // -------------------------------------------------------------
55 
56 
57 // Calculate and invert the diagonal elements once and for all
58 // -------------------------------------------------------------
59  valarray<double> invDIAG(A.dim1());
60  for (int i(0); i < invDIAG.size(); ++i) invDIAG[i] = 1.0 / A(i,i);
61 
62  valarray<complex<double> > xold(xk);
63  int iteration(0); // used to count iterations
64  int conv(0); // used to test convergence
65 
66 // Start the iteration loop
67 // -------------------------------------------------------------
68  while ( (iteration++ < MAXiter) && (conv < b.size()) ) {
69 
70  xold = xk;
71  for (int i(0); i < A.dim1(); ++i){
72  complex<double> sigma(0.0); // Temporary sum
73  for (int j(0); j < i; ++j){
74  sigma += A(i,j)*xk[j];
75  }
76  for (int j(i+1); j < A.dim2(); ++j){
77  sigma += A(i,j)*xk[j];
78  }
79  xk[i] = invDIAG[i] * (b[i] - sigma);
80  }
81 
82  // Calculate Dx = x_old - x_new
83  xold -= xk;
84 
85 // If the relative error < prescribed tolerance everywhere the method has converged
86 // |Dx(i)| < t*|x(i)| + eps
87  conv = 0;
88  while ( ( conv < b.size() ) &&
89  ( abs(xold[conv]) < (tol*abs(xk[conv] + 20.0*DBL_MIN)) ) ){
90  ++conv;
91  }
92 
93  //----> Output for testing
94  //--------------------------------
95  //cout << "iteration = " << iteration << " ";
96  //for (int i(0); i < b.size(); ++i){
97  // cout << xk[i] << " ";
98  //}
99  //cout << "\n";
100  //--------------------------------
101 
102  }
103 
104  //----> Output for testing
105  //--------------------------------
106  // cout << "Iterations = " << iteration-1 <<"\n";
107  //for (int i(0); i < b.size(); ++i) {
108  // cout << "Error |Dx| = " << abs(xold[i])
109  // << ", "
110  // << "Tolerance * |x| = " << tol*abs(xk[i]) <<"\n";
111  //}
112  //--------------------------------
113 
114  return true;
115  }
Array2D::dim1
size_t dim1() const
Definition: lib-array.h:289
Array2D::dim2
size_t dim2() const
Definition: lib-array.h:290
Here is the call graph for this function:
Here is the caller graph for this function:

◆ Thomas_Tridiagonal() [1/3]

bool Thomas_Tridiagonal ( Array2D< double > &  A,
valarray< double > &  d,
valarray< double > &  xk 
)

The tridiagonal solver for implicit collisions.

Parameters
AInput matrix
dRight Side
xkSolution vector
Returns
xk is returned as the solution vector

Definition at line 216 of file nmethods.cpp.

References Array2D< T >::dim1(), Array2D< T >::dim2(), and TridiagonalSolve().

Referenced by interspecies_flm_implicit_step::advance(), self_flm_implicit_step::advance(), Magnetic_Field_1D_f1::implicit(), Magnetic_Field_1D::implicit(), self_f00_implicit_collisions::loop(), and self_f00_implicit_step::takestep().

218  {
219 //-------------------------------------------------------------------
220 // Fills solution into xk. The other matrices are not modified
221 // The function returns "false" if the matrix A is not diagonally
222 // dominant
223 //-------------------------------------------------------------------
224 
225 // The Matrices all have the right dimensions
226 // -------------------------------------------------------------
227  if ( ( A.dim1() != A.dim2() ) ||
228  ( A.dim1() != d.size() ) ||
229  ( A.dim1() != xk.size() ) ) {
230  cout << "Error: The Matrices don't have the right dimensions!" << endl;
231  exit(1);
232  }
233 // -------------------------------------------------------------
234 
235  valarray<double> a(d.size()), b(d.size()), c(d.size());
236 
237  for (int i(0); i < A.dim1()-1; ++i){
238  a[i+1] = A(i+1,i);
239  }
240  for (int i(0); i < A.dim1(); ++i){
241  b[i] = A(i,i);
242  }
243  for (int i(0); i < A.dim1()-1; ++i){
244  c[i] = A(i,i+1);
245  }
246 
247 // valarray< double > dcopy(d);
248  TridiagonalSolve(a,b,c,d,xk);
249 
250  return true;
251 }
Array2D::dim1
size_t dim1() const
Definition: lib-array.h:289
Array2D::dim2
size_t dim2() const
Definition: lib-array.h:290
TridiagonalSolve
void TridiagonalSolve(const valarray< double > &a, const valarray< double > &b, valarray< double > &c, valarray< complex< double > > d, valarray< complex< double > > &x)
Set up tridiagonal solver.
Definition: nmethods.cpp:121
Here is the call graph for this function:
Here is the caller graph for this function:

◆ Thomas_Tridiagonal() [2/3]

bool Thomas_Tridiagonal ( Array2D< double > &  A,
valarray< complex< double > > &  d,
valarray< complex< double > > &  xk 
)

Definition at line 254 of file nmethods.cpp.

References Array2D< T >::dim1(), Array2D< T >::dim2(), and TridiagonalSolve().

256  {
257 //-------------------------------------------------------------------
258 // Fills solution into xk. The other matrices are not modified
259 // The function returns "false" if the matrix A is not diagonally
260 // dominant
261 //-------------------------------------------------------------------
262 
263 // The Matrices all have the right dimensions
264 // -------------------------------------------------------------
265  if ( ( A.dim1() != A.dim2() ) ||
266  ( A.dim1() != d.size() ) ||
267  ( A.dim1() != xk.size() ) ) {
268  cout << "Error: The Matrices don't have the right dimensions!" << endl;
269  exit(1);
270  }
271 // -------------------------------------------------------------
272 
273  valarray<double> a(d.size()), b(d.size()), c(d.size());
274 
275  for (int i(0); i < A.dim1()-1; ++i){
276  a[i+1] = A(i+1,i);
277  }
278  for (int i(0); i < A.dim1(); ++i){
279  b[i] = A(i,i);
280  }
281  for (int i(0); i < A.dim1()-1; ++i){
282  c[i] = A(i,i+1);
283  }
284 
285 // valarray< double > dcopy(d);
286  TridiagonalSolve(a,b,c,d,xk);
287 
288  return true;
289 }
Array2D::dim1
size_t dim1() const
Definition: lib-array.h:289
Array2D::dim2
size_t dim2() const
Definition: lib-array.h:290
TridiagonalSolve
void TridiagonalSolve(const valarray< double > &a, const valarray< double > &b, valarray< double > &c, valarray< complex< double > > d, valarray< complex< double > > &x)
Set up tridiagonal solver.
Definition: nmethods.cpp:121
Here is the call graph for this function:

◆ Thomas_Tridiagonal() [3/3]

bool Thomas_Tridiagonal ( Array2D< complex< double > > &  A,
valarray< complex< double > > &  d,
valarray< complex< double > > &  xk 
)

◆ TridiagonalSolve() [1/3]

void TridiagonalSolve ( const valarray< double > &  a,
const valarray< double > &  b,
valarray< double > &  c,
valarray< complex< double > >  d,
valarray< complex< double > > &  x 
)

Set up tridiagonal solver.

Parameters
[in]a???
[in]b???
c???
[in]dright side
xsolution

Definition at line 121 of file nmethods.cpp.

Referenced by Thomas_Tridiagonal().

125  {
126 //-------------------------------------------------------------------
127 // Fills solution into x. Warning: will modify c and d!
128 //-------------------------------------------------------------------
129  size_t n(x.size());
130  // Modify the coefficients.
131  c[0] /= b[0]; // Division by zero risk.
132  d[0] /= b[0]; // Division by zero would imply a singular matrix.
133  for (int i(1); i < n; ++i){
134  double id(1.0/(b[i]-c[i-1]*a[i])); // Division by zero risk.
135  c[i] *= id; // Last value calculated is redundant.
136  d[i] -= d[i-1] * a[i];
137  d[i] *= id; // d[i] = (d[i] - d[i-1] * a[i]) * id
138  }
139 
140  // Now back substitute.
141  x[n-1] = d[n-1];
142  for (int i(n-2); i > -1; --i){
143  x[i] = d[i];
144  x[i] -= c[i] * x[i+1]; // x[i] = d[i] - c[i] * x[i + 1];
145  }
146  }
Here is the caller graph for this function:

◆ TridiagonalSolve() [2/3]

void TridiagonalSolve ( const valarray< double > &  a,
const valarray< double > &  b,
valarray< double > &  c,
valarray< double >  d,
valarray< double > &  x 
)

Definition at line 151 of file nmethods.cpp.

155  {
156 //-------------------------------------------------------------------
157 // Fills solution into x. Warning: will modify c and d!
158 //-------------------------------------------------------------------
159  size_t n(x.size());
160  // Modify the coefficients.
161  c[0] /= b[0]; // Division by zero risk.
162  d[0] /= b[0]; // Division by zero would imply a singular matrix.
163  for (int i(1); i < n; ++i){
164  double id(1.0/(b[i]-c[i-1]*a[i])); // Division by zero risk.
165  c[i] *= id; // Last value calculated is redundant.
166  d[i] -= d[i-1] * a[i];
167  d[i] *= id; // d[i] = (d[i] - d[i-1] * a[i]) * id
168  }
169 
170  // Now back substitute.
171  x[n-1] = d[n-1];
172  for (int i(n-2); i > -1; --i){
173  x[i] = d[i];
174  x[i] -= c[i] * x[i+1]; // x[i] = d[i] - c[i] * x[i + 1];
175  }
176 }

◆ TridiagonalSolve() [3/3]

void TridiagonalSolve ( const valarray< complex< double > > &  a,
const valarray< complex< double > > &  b,
valarray< complex< double > > &  c,
valarray< complex< double > >  d,
valarray< complex< double > > &  x 
)

◆ vfloat() [1/3]

valarray<float> vfloat ( const valarray< double > &  vDouble)

Definition at line 381 of file nmethods.cpp.

Referenced by Output_Data::Output_Preprocessor_1D::Jx(), Output_Data::Output_Preprocessor_1D::Jy(), Output_Data::Output_Preprocessor_1D::n(), Output_Data::Output_Preprocessor_1D::Qx(), Output_Data::Output_Preprocessor_1D::Qy(), Output_Data::Output_Preprocessor_1D::T(), Output_Data::Output_Preprocessor_1D::vNx(), Output_Data::Output_Preprocessor_1D::vNy(), and Output_Data::Output_Preprocessor_1D::vNz().

381  {
382  valarray<float> vf(vDouble.size());
383  for (size_t i(0); i < vf.size(); ++i) {
384  vf[i] = static_cast<float>(vDouble[i]);
385  }
386  return vf;
387 }
Here is the caller graph for this function:

◆ vfloat() [2/3]

vector<float> vfloat ( const vector< double >  vDouble)

Convert data structure to float structure.

Parameters
[in]vDoubleThe v double
Returns
float structure

Definition at line 397 of file nmethods.cpp.

397  {
398  vector<float> vf;
399  for (size_t i(0); i < vDouble.size(); ++i) {
400  vf.push_back(static_cast<float>(vDouble[i]));
401  }
402  return vf;
403 }

◆ vfloat() [3/3]

vector<float> vfloat ( const vector< complex< double > >  vDouble)

Definition at line 404 of file nmethods.cpp.

404  {
405  vector<float> vf;
406  for (size_t i(0); i < vDouble.size(); ++i) {
407  vf.push_back(static_cast<float>(vDouble[i].real()));
408  }
409  return vf;
410 }

◆ vfloat_complex()

vector<float> vfloat_complex ( const vector< complex< double > >  vDouble)

Definition at line 412 of file nmethods.cpp.

Referenced by Output_Data::Output_Preprocessor_1D::Jz(), Output_Data::Output_Preprocessor_1D::Qz(), and Output_Data::Output_Preprocessor_1D::vNz().

412  {
413  vector<float> vf;
414  for (size_t i(0); i < vDouble.size(); ++i) {
415  vf.push_back(static_cast<float>(vDouble[i].imag()));
416  }
417  return vf;
418 }
Here is the caller graph for this function: