From 7fa9c70b863eba9d4c37393fbf1d60b3f4a9be72 Mon Sep 17 00:00:00 2001 From: John Date: Wed, 2 Mar 2022 17:14:08 -0500 Subject: [PATCH 1/7] ported field perturbation still need particles --- src/psc_harris_yz.cxx | 42 ++++++++++++++++++++++++++++++++++-------- 1 file changed, 34 insertions(+), 8 deletions(-) diff --git a/src/psc_harris_yz.cxx b/src/psc_harris_yz.cxx index 940820c267..213a371126 100644 --- a/src/psc_harris_yz.cxx +++ b/src/psc_harris_yz.cxx @@ -49,6 +49,7 @@ struct PscHarrisParams double background_Ti; double bg; // guide field as fraction of B0 + double dbz_b0; // The following parameters are calculated from the above / and other // information @@ -57,6 +58,14 @@ struct PscHarrisParams double b0; // B0 double L; // sheet width in d_e double TTi, TTe; + double wpe_wce; + + double wpe, wpi, wce, wci; + double Lx, Ly, Lz; // size of box + double Lpert_Lx; + double Lpert; // wavelength of perturbation + double dbx; // Perturbation in Bz relative to Bo (Only change here) + double dbz; // Set Bx perturbation so that div(B) = 0 }; // ====================================================================== @@ -146,8 +155,27 @@ void setupParameters() g.background_n = .002; g.background_Te = .001; g.background_Ti = .001; + g.wpe_wce = 2.; g.bg = 0.; + g.dbz_b0 = .03; + g.Lpert_Lx = 1.; + double c = 1.; // Speed of light + g.wci = 1. / (g.mass_ratio * g.wpe_wce); // Ion cyclotron frequency + g.wce = g.wci * g.mass_ratio; // Electron cyclotron freqeuncy + g.wpe = g.wce * g.wpe_wce; // electron plasma frequency + + g.wpi = g.wpe / sqrt(g.mass_ratio); // ion plasma frequency + g.d_i = c / g.wpi; // ion inertial length + g.Lx = g.Lx_di * g.d_i; // size of box in x dimension + g.Ly = g.Ly_di * g.d_i; // size of box in y dimension + g.Lz = g.Lz_di * g.d_i; // size of box in z dimension + + g.Lpert = g.Lpert_Lx * g.Lx; // wavelength of perturbation + g.dbz = + g.dbz_b0 * g.b0; // Perturbation in Bz relative to Bo (Only change here) + g.dbx = + -g.dbz * g.Lpert / (2. * g.Lz); // Set Bx perturbation so that div(B) = 0 } // ====================================================================== @@ -276,8 +304,8 @@ void initializeParticles(SetupParticles& setup_particles, void initializeFields(MfieldsState& mflds) { - double b0 = g.b0; //, dbx = g.dbx, dbz = g.dbz; - double L = g.L; //, Ly = g.Ly, Lz = g.Lz;//, Lpert = g.Lpert; + double b0 = g.b0, dbx = g.dbx, dbz = g.dbz; + double L = g.L, Lx = g.Lx, Ly = g.Ly, Lz = g.Lz, Lpert = g.Lpert; // double cs = cos(g.theta), sn = sin(g.theta); setupFields(mflds, [&](int m, double crd[3]) { @@ -285,16 +313,14 @@ void initializeFields(MfieldsState& mflds) switch (m) { case HY: - return b0 * - tanh(z / L); // + - // dbx * cos(2. * M_PI * (x - .5 * Lx) / Lpert) * - // sin(M_PI * z / Lz); + return b0 * tanh(z / L); + +dbx* cos(2. * M_PI * (x - .5 * Lx) / Lpert) * sin(M_PI * z / Lz); case HX: return b0 * g.bg; case HZ: - // return dbz * cos(M_PI * z / Lz) * - // sin(2.0 * M_PI * (x - 0.5 * Lx) / Lpert); + return dbz * cos(M_PI * z / Lz) * + sin(2.0 * M_PI * (x - 0.5 * Lx) / Lpert); return 0.; case JYI: return 0.; // FIXME From 8818688590daa87dcbf37d81eb3988f24ec76716 Mon Sep 17 00:00:00 2001 From: John Date: Fri, 4 Mar 2022 10:52:23 -0500 Subject: [PATCH 2/7] fixed initial field conditions --- src/psc_harris_yz.cxx | 72 +++++++++++++++++++++++-------------------- 1 file changed, 38 insertions(+), 34 deletions(-) diff --git a/src/psc_harris_yz.cxx b/src/psc_harris_yz.cxx index 213a371126..6ffe06dc82 100644 --- a/src/psc_harris_yz.cxx +++ b/src/psc_harris_yz.cxx @@ -44,27 +44,29 @@ struct PscHarrisParams double mass_ratio; double lambda0; - double background_n; - double background_Te; - double background_Ti; + //double background_n; + //double background_Te; + //double background_Ti; double bg; // guide field as fraction of B0 + double theta; double dbz_b0; + double Lpert_Ly; // The following parameters are calculated from the above / and other // information double d_i; double b0; // B0 + //double n0; double L; // sheet width in d_e double TTi, TTe; double wpe_wce; double wpe, wpi, wce, wci; double Lx, Ly, Lz; // size of box - double Lpert_Lx; double Lpert; // wavelength of perturbation - double dbx; // Perturbation in Bz relative to Bo (Only change here) + double dby; // Perturbation in Bz relative to Bo (Only change here) double dbz; // Set Bx perturbation so that div(B) = 0 }; @@ -141,40 +143,48 @@ void setupParameters() // read_checkpoint_filename = "checkpoint_500.bp"; // -- Set some parameters specific to this case + + // PIC units, only used in this scope + //double eps0 = 1; + double me = 1; + double ec = 1; + double c = 1; // Speed of light + g.Lx_di = 1.; g.Ly_di = 40.; g.Lz_di = 10.; - g.L_di = .5; + g.Lpert_Ly = 1.; g.BB = 0.; g.Zi = 1.; g.mass_ratio = 100.; g.lambda0 = 20.; + g.bg = 0.; + g.theta = 0; + g.dbz_b0 = .03; - g.background_n = .002; - g.background_Te = .001; - g.background_Ti = .001; + g.TTe = .2; + g.TTi = .1; g.wpe_wce = 2.; - g.bg = 0.; - g.dbz_b0 = .03; - g.Lpert_Lx = 1.; - double c = 1.; // Speed of light g.wci = 1. / (g.mass_ratio * g.wpe_wce); // Ion cyclotron frequency g.wce = g.wci * g.mass_ratio; // Electron cyclotron freqeuncy g.wpe = g.wce * g.wpe_wce; // electron plasma frequency - g.wpi = g.wpe / sqrt(g.mass_ratio); // ion plasma frequency + g.d_i = c / g.wpi; // ion inertial length + g.L = g.L_di * g.d_i; g.Lx = g.Lx_di * g.d_i; // size of box in x dimension g.Ly = g.Ly_di * g.d_i; // size of box in y dimension g.Lz = g.Lz_di * g.d_i; // size of box in z dimension - g.Lpert = g.Lpert_Lx * g.Lx; // wavelength of perturbation + g.b0 = me * c * g.wce / ec; // Asymptotic magnetic field strength + //g.n0 = me * eps0 * wpe * wpe / (ec * ec); // Peak electron (ion) density + g.Lpert = g.Lpert_Ly * g.Ly; // wavelength of perturbation g.dbz = g.dbz_b0 * g.b0; // Perturbation in Bz relative to Bo (Only change here) - g.dbx = + g.dby = -g.dbz * g.Lpert / (2. * g.Lz); // Set Bx perturbation so that div(B) = 0 } @@ -193,11 +203,9 @@ Grid_t* setupGrid() kinds[MY_ION] = {g.Zi, g.mass_ratio * g.Zi, "i"}; kinds[MY_ELECTRON] = {-1., 1., "e"}; - g.d_i = sqrt(kinds[MY_ION].m / kinds[MY_ION].q); - mpi_printf(MPI_COMM_WORLD, "d_e = %g, d_i = %g\n", 1., g.d_i); mpi_printf(MPI_COMM_WORLD, "lambda_De (background) = %g\n", - sqrt(g.background_Te)); + sqrt(g.TTe)); // --- setup domain Grid_t::Real3 LL = {g.Lx_di * g.d_i, g.Ly_di * g.d_i, @@ -221,11 +229,7 @@ Grid_t* setupGrid() mpi_printf(MPI_COMM_WORLD, "dt = %g\n", dt); - g.b0 = 1.; // FIXME - g.L = g.L_di * g.d_i; - g.TTe = .2; - g.TTi = .1; Int3 ibn = {2, 2, 2}; if (Dim::InvarX::value) { @@ -304,26 +308,26 @@ void initializeParticles(SetupParticles& setup_particles, void initializeFields(MfieldsState& mflds) { - double b0 = g.b0, dbx = g.dbx, dbz = g.dbz; - double L = g.L, Lx = g.Lx, Ly = g.Ly, Lz = g.Lz, Lpert = g.Lpert; - // double cs = cos(g.theta), sn = sin(g.theta); + double b0 = g.b0, dby = g.dby, dbz = g.dbz; + double L = g.L, Ly = g.Ly, Lz = g.Lz, Lpert = g.Lpert; + double cs = cos(g.theta), sn = sin(g.theta); setupFields(mflds, [&](int m, double crd[3]) { - double x = crd[1], z = crd[2]; + double y = crd[1], z = crd[2]; switch (m) { + case HX: return -sn * b0 * tanh(z / L) + b0 * g.bg; + case HY: - return b0 * tanh(z / L); - +dbx* cos(2. * M_PI * (x - .5 * Lx) / Lpert) * sin(M_PI * z / Lz); - - case HX: return b0 * g.bg; + return cs * b0 * tanh(z / L) + + dby * cos(2. * M_PI * (y - .5 * Ly) / Lpert) * + sin(M_PI * z / Lz); case HZ: return dbz * cos(M_PI * z / Lz) * - sin(2.0 * M_PI * (x - 0.5 * Lx) / Lpert); - return 0.; + sin(2.0 * M_PI * (y - 0.5 * Ly) / Lpert); - case JYI: return 0.; // FIXME + case JXI: return 0.; // FIXME default: return 0.; } From e58a6a5ea9c95a8eeb4e88a7240d3b146f89af53 Mon Sep 17 00:00:00 2001 From: John Date: Fri, 4 Mar 2022 11:00:55 -0500 Subject: [PATCH 3/7] calculated prt TTe TTi --- src/psc_harris_yz.cxx | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/psc_harris_yz.cxx b/src/psc_harris_yz.cxx index 6ffe06dc82..c738231ea1 100644 --- a/src/psc_harris_yz.cxx +++ b/src/psc_harris_yz.cxx @@ -42,6 +42,7 @@ struct PscHarrisParams double BB; double Zi; double mass_ratio; + double Ti_Te; double lambda0; //double background_n; @@ -149,6 +150,7 @@ void setupParameters() double me = 1; double ec = 1; double c = 1; // Speed of light + double eps0 = 1; g.Lx_di = 1.; g.Ly_di = 40.; @@ -159,13 +161,14 @@ void setupParameters() g.BB = 0.; g.Zi = 1.; g.mass_ratio = 100.; + g.Ti_Te = 5.; g.lambda0 = 20.; g.bg = 0.; g.theta = 0; g.dbz_b0 = .03; - g.TTe = .2; - g.TTi = .1; + g.TTe = me * sqr(c) / (2. * eps0 * sqr(g.wpe_wce) * (1. + g.Ti_Te)); + g.TTi = g.TTe * g.Ti_Te; g.wpe_wce = 2.; g.wci = 1. / (g.mass_ratio * g.wpe_wce); // Ion cyclotron frequency From 18deac43f0354b4f936124aca118ed94b34e991f Mon Sep 17 00:00:00 2001 From: John Date: Fri, 4 Mar 2022 14:49:37 -0500 Subject: [PATCH 4/7] same params as xz for comparison, initial conditions look good, evolution is still off --- src/psc_harris_yz.cxx | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/src/psc_harris_yz.cxx b/src/psc_harris_yz.cxx index c738231ea1..e7b5d361d8 100644 --- a/src/psc_harris_yz.cxx +++ b/src/psc_harris_yz.cxx @@ -11,7 +11,7 @@ #ifdef USE_CUDA #include "cuda_bits.h" #endif - +#define PRINT(c) {std::cout << "JOHN " << c << std::endl;} // ====================================================================== // Particle kinds // @@ -154,22 +154,22 @@ void setupParameters() g.Lx_di = 1.; g.Ly_di = 40.; - g.Lz_di = 10.; + g.Lz_di = 20.; g.L_di = .5; g.Lpert_Ly = 1.; g.BB = 0.; g.Zi = 1.; - g.mass_ratio = 100.; + g.mass_ratio = 16.; g.Ti_Te = 5.; g.lambda0 = 20.; g.bg = 0.; g.theta = 0; g.dbz_b0 = .03; + g.wpe_wce = 2.; g.TTe = me * sqr(c) / (2. * eps0 * sqr(g.wpe_wce) * (1. + g.Ti_Te)); g.TTi = g.TTe * g.Ti_Te; - g.wpe_wce = 2.; g.wci = 1. / (g.mass_ratio * g.wpe_wce); // Ion cyclotron frequency g.wce = g.wci * g.mass_ratio; // Electron cyclotron freqeuncy @@ -213,8 +213,10 @@ Grid_t* setupGrid() // --- setup domain Grid_t::Real3 LL = {g.Lx_di * g.d_i, g.Ly_di * g.d_i, g.Lz_di * g.d_i}; // domain size (in d_e) - Int3 gdims = {1, 512, 128}; - Int3 np = {1, 4, 1}; + //Int3 gdims = {1, 512, 128}; + //Int3 np = {1, 4, 1}; + Int3 gdims = {1, 640, 320}; + Int3 np = {1, 20, 10}; Grid_t::Domain domain{gdims, LL, -.5 * LL, np}; @@ -225,7 +227,7 @@ Grid_t* setupGrid() // -- setup normalization auto norm_params = Grid_t::NormalizationParams::dimensionless(); - norm_params.nicell = 100; + norm_params.nicell = 10000; double dt = psc_params.cfl * courant_length(domain); Grid_t::Normalization norm{norm_params}; From f79354200d3bf7394fde0be2af62089f915efa84 Mon Sep 17 00:00:00 2001 From: John Date: Fri, 4 Mar 2022 14:50:33 -0500 Subject: [PATCH 5/7] clang format --- src/psc_harris_yz.cxx | 54 +++++++++++++++++++++---------------------- 1 file changed, 27 insertions(+), 27 deletions(-) diff --git a/src/psc_harris_yz.cxx b/src/psc_harris_yz.cxx index e7b5d361d8..49c3ce5bad 100644 --- a/src/psc_harris_yz.cxx +++ b/src/psc_harris_yz.cxx @@ -11,7 +11,10 @@ #ifdef USE_CUDA #include "cuda_bits.h" #endif -#define PRINT(c) {std::cout << "JOHN " << c << std::endl;} +#define PRINT(c) \ + { \ + std::cout << "JOHN " << c << std::endl; \ + } // ====================================================================== // Particle kinds // @@ -45,9 +48,9 @@ struct PscHarrisParams double Ti_Te; double lambda0; - //double background_n; - //double background_Te; - //double background_Ti; + // double background_n; + // double background_Te; + // double background_Ti; double bg; // guide field as fraction of B0 double theta; @@ -59,16 +62,16 @@ struct PscHarrisParams double d_i; double b0; // B0 - //double n0; - double L; // sheet width in d_e + // double n0; + double L; // sheet width in d_e double TTi, TTe; double wpe_wce; double wpe, wpi, wce, wci; double Lx, Ly, Lz; // size of box - double Lpert; // wavelength of perturbation - double dby; // Perturbation in Bz relative to Bo (Only change here) - double dbz; // Set Bx perturbation so that div(B) = 0 + double Lpert; // wavelength of perturbation + double dby; // Perturbation in Bz relative to Bo (Only change here) + double dbz; // Set Bx perturbation so that div(B) = 0 }; // ====================================================================== @@ -144,12 +147,12 @@ void setupParameters() // read_checkpoint_filename = "checkpoint_500.bp"; // -- Set some parameters specific to this case - + // PIC units, only used in this scope - //double eps0 = 1; + // double eps0 = 1; double me = 1; double ec = 1; - double c = 1; // Speed of light + double c = 1; // Speed of light double eps0 = 1; g.Lx_di = 1.; @@ -174,16 +177,16 @@ void setupParameters() g.wci = 1. / (g.mass_ratio * g.wpe_wce); // Ion cyclotron frequency g.wce = g.wci * g.mass_ratio; // Electron cyclotron freqeuncy g.wpe = g.wce * g.wpe_wce; // electron plasma frequency - g.wpi = g.wpe / sqrt(g.mass_ratio); // ion plasma frequency + g.wpi = g.wpe / sqrt(g.mass_ratio); // ion plasma frequency - g.d_i = c / g.wpi; // ion inertial length + g.d_i = c / g.wpi; // ion inertial length g.L = g.L_di * g.d_i; - g.Lx = g.Lx_di * g.d_i; // size of box in x dimension - g.Ly = g.Ly_di * g.d_i; // size of box in y dimension - g.Lz = g.Lz_di * g.d_i; // size of box in z dimension + g.Lx = g.Lx_di * g.d_i; // size of box in x dimension + g.Ly = g.Ly_di * g.d_i; // size of box in y dimension + g.Lz = g.Lz_di * g.d_i; // size of box in z dimension - g.b0 = me * c * g.wce / ec; // Asymptotic magnetic field strength - //g.n0 = me * eps0 * wpe * wpe / (ec * ec); // Peak electron (ion) density + g.b0 = me * c * g.wce / ec; // Asymptotic magnetic field strength + // g.n0 = me * eps0 * wpe * wpe / (ec * ec); // Peak electron (ion) density g.Lpert = g.Lpert_Ly * g.Ly; // wavelength of perturbation g.dbz = g.dbz_b0 * g.b0; // Perturbation in Bz relative to Bo (Only change here) @@ -207,14 +210,13 @@ Grid_t* setupGrid() kinds[MY_ELECTRON] = {-1., 1., "e"}; mpi_printf(MPI_COMM_WORLD, "d_e = %g, d_i = %g\n", 1., g.d_i); - mpi_printf(MPI_COMM_WORLD, "lambda_De (background) = %g\n", - sqrt(g.TTe)); + mpi_printf(MPI_COMM_WORLD, "lambda_De (background) = %g\n", sqrt(g.TTe)); // --- setup domain Grid_t::Real3 LL = {g.Lx_di * g.d_i, g.Ly_di * g.d_i, g.Lz_di * g.d_i}; // domain size (in d_e) - //Int3 gdims = {1, 512, 128}; - //Int3 np = {1, 4, 1}; + // Int3 gdims = {1, 512, 128}; + // Int3 np = {1, 4, 1}; Int3 gdims = {1, 640, 320}; Int3 np = {1, 20, 10}; @@ -234,8 +236,6 @@ Grid_t* setupGrid() mpi_printf(MPI_COMM_WORLD, "dt = %g\n", dt); - - Int3 ibn = {2, 2, 2}; if (Dim::InvarX::value) { ibn[0] = 0; @@ -322,10 +322,10 @@ void initializeFields(MfieldsState& mflds) switch (m) { case HX: return -sn * b0 * tanh(z / L) + b0 * g.bg; - + case HY: return cs * b0 * tanh(z / L) + - dby * cos(2. * M_PI * (y - .5 * Ly) / Lpert) * + dby * cos(2. * M_PI * (y - .5 * Ly) / Lpert) * sin(M_PI * z / Lz); case HZ: From 530f92655794ce959782df7a32711df27320e4fe Mon Sep 17 00:00:00 2001 From: John Date: Mon, 7 Mar 2022 10:38:04 -0500 Subject: [PATCH 6/7] added background species and conducting wall/prt reflecting boundary conditions --- src/psc_harris_yz.cxx | 64 ++++++++++++++++++++----------------------- 1 file changed, 30 insertions(+), 34 deletions(-) diff --git a/src/psc_harris_yz.cxx b/src/psc_harris_yz.cxx index 49c3ce5bad..c718ba52bc 100644 --- a/src/psc_harris_yz.cxx +++ b/src/psc_harris_yz.cxx @@ -30,6 +30,8 @@ enum { MY_ELECTRON, MY_ION, + MY_ELECTRON_BG, + MY_ION_BG, N_MY_KINDS, }; @@ -46,8 +48,9 @@ struct PscHarrisParams double Zi; double mass_ratio; double Ti_Te; + double Tib_Ti, Teb_Te; double lambda0; - + double nb_n0; //background density // double background_n; // double background_Te; // double background_Ti; @@ -165,7 +168,10 @@ void setupParameters() g.Zi = 1.; g.mass_ratio = 16.; g.Ti_Te = 5.; + g.Tib_Ti = 0.333; + g.Teb_Te = 0.333; g.lambda0 = 20.; + g.nb_n0 = 0.05; g.bg = 0.; g.theta = 0; g.dbz_b0 = .03; @@ -208,6 +214,8 @@ Grid_t* setupGrid() Grid_t::Kinds kinds(N_MY_KINDS); kinds[MY_ION] = {g.Zi, g.mass_ratio * g.Zi, "i"}; kinds[MY_ELECTRON] = {-1., 1., "e"}; + kinds[MY_ION_BG] = {g.Zi, g.mass_ratio * g.Zi, "i_bg"}; + kinds[MY_ELECTRON_BG] = {-1., 1., "e_bg"}; mpi_printf(MPI_COMM_WORLD, "d_e = %g, d_i = %g\n", 1., g.d_i); mpi_printf(MPI_COMM_WORLD, "lambda_De (background) = %g\n", sqrt(g.TTe)); @@ -215,17 +223,16 @@ Grid_t* setupGrid() // --- setup domain Grid_t::Real3 LL = {g.Lx_di * g.d_i, g.Ly_di * g.d_i, g.Lz_di * g.d_i}; // domain size (in d_e) - // Int3 gdims = {1, 512, 128}; - // Int3 np = {1, 4, 1}; + Int3 gdims = {1, 640, 320}; Int3 np = {1, 20, 10}; Grid_t::Domain domain{gdims, LL, -.5 * LL, np}; - psc::grid::BC bc{{BND_FLD_PERIODIC, BND_FLD_PERIODIC, BND_FLD_PERIODIC}, - {BND_FLD_PERIODIC, BND_FLD_PERIODIC, BND_FLD_PERIODIC}, - {BND_PRT_PERIODIC, BND_PRT_PERIODIC, BND_PRT_PERIODIC}, - {BND_PRT_PERIODIC, BND_PRT_PERIODIC, BND_PRT_PERIODIC}}; + psc::grid::BC bc{{BND_FLD_PERIODIC, BND_FLD_PERIODIC, BND_FLD_CONDUCTING_WALL}, + {BND_FLD_PERIODIC, BND_FLD_PERIODIC, BND_FLD_CONDUCTING_WALL}, + {BND_PRT_PERIODIC, BND_PRT_PERIODIC, BND_PRT_REFLECTING}, + {BND_PRT_PERIODIC, BND_PRT_PERIODIC, BND_PRT_REFLECTING}}; // -- setup normalization auto norm_params = Grid_t::NormalizationParams::dimensionless(); @@ -259,15 +266,6 @@ void initializeParticles(SetupParticles& setup_particles, // -- set particle initial condition partitionAndSetupParticles(setup_particles, balance, grid_ptr, mprts, [&](int kind, Double3 crd, psc_particle_npt& npt) { - // struct psc_harris* harris = - // to_psc_harris(psc); - - // double BB = harris->B0; - // double LLL = harris->LLL; - // double nnb = harris->nb; - // double TTi = harris->Ti; - // double TTe = harris->Te; - switch (kind) { case MY_ION: // ion drifting npt.n = 1. / sqr(cosh(crd[2] / g.L)); @@ -277,15 +275,14 @@ void initializeParticles(SetupParticles& setup_particles, npt.T[2] = g.TTi; npt.kind = MY_ION; break; - // case 1: // ion bg - // npt->n = nnb; - // npt->q = 1.; - // npt->m = harris->mi_over_me; - // npt->T[0] = TTi; - // npt->T[1] = TTi; - // npt->T[2] = TTi; - // npt->kind = KIND_ION; - // break; + case MY_ION_BG: // ion bg + npt.n = g.nb_n0 / sqr(cosh(crd[2] / g.L)); + npt.p[0] = 2. * g.Tib_Ti * g.TTi / g.b0 / g.L; + npt.T[0] = g.Tib_Ti * g.TTi; + npt.T[1] = g.Tib_Ti * g.TTi; + npt.T[2] = g.Tib_Ti * g.TTi; + npt.kind = MY_ION_BG; + break; case MY_ELECTRON: // electron drifting npt.n = 1. / sqr(cosh(crd[2] / g.L)); npt.p[0] = -2. * g.TTe / g.b0 / g.L; @@ -294,15 +291,14 @@ void initializeParticles(SetupParticles& setup_particles, npt.T[2] = g.TTe; npt.kind = MY_ELECTRON; break; - // case 3: // electron bg - // npt->n = nnb; - // npt->q = -1.; - // npt->m = 1.; - // npt->T[0] = TTe; - // npt->T[1] = TTe; - // npt->T[2] = TTe; - // npt->kind = KIND_ELECTRON; - // break; + case MY_ELECTRON_BG: // electron bg + npt.n = g.nb_n0 / sqr(cosh(crd[2] / g.L)); + npt.p[0] = -2. * g.Teb_Te * g.TTe / g.b0 / g.L; + npt.T[0] = g.Teb_Te * g.TTe; + npt.T[1] = g.Teb_Te * g.TTe; + npt.T[2] = g.Teb_Te * g.TTe; + npt.kind = MY_ELECTRON_BG; + break; default: assert(0); } }); From b76dd66cda17092d1ff9def5c5bfdee10a691a05 Mon Sep 17 00:00:00 2001 From: John Date: Mon, 7 Mar 2022 15:09:03 -0500 Subject: [PATCH 7/7] corrected drift velocity sign --- src/psc_harris_yz.cxx | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/psc_harris_yz.cxx b/src/psc_harris_yz.cxx index c718ba52bc..8250a65c8e 100644 --- a/src/psc_harris_yz.cxx +++ b/src/psc_harris_yz.cxx @@ -269,7 +269,7 @@ void initializeParticles(SetupParticles& setup_particles, switch (kind) { case MY_ION: // ion drifting npt.n = 1. / sqr(cosh(crd[2] / g.L)); - npt.p[0] = 2. * g.TTi / g.b0 / g.L; + npt.p[0] = -2. * g.TTi / g.b0 / g.L; npt.T[0] = g.TTi; npt.T[1] = g.TTi; npt.T[2] = g.TTi; @@ -277,7 +277,7 @@ void initializeParticles(SetupParticles& setup_particles, break; case MY_ION_BG: // ion bg npt.n = g.nb_n0 / sqr(cosh(crd[2] / g.L)); - npt.p[0] = 2. * g.Tib_Ti * g.TTi / g.b0 / g.L; + npt.p[0] = 0.;//2. * g.Tib_Ti * g.TTi / g.b0 / g.L; npt.T[0] = g.Tib_Ti * g.TTi; npt.T[1] = g.Tib_Ti * g.TTi; npt.T[2] = g.Tib_Ti * g.TTi; @@ -285,7 +285,7 @@ void initializeParticles(SetupParticles& setup_particles, break; case MY_ELECTRON: // electron drifting npt.n = 1. / sqr(cosh(crd[2] / g.L)); - npt.p[0] = -2. * g.TTe / g.b0 / g.L; + npt.p[0] = 2. * g.TTe / g.b0 / g.L; npt.T[0] = g.TTe; npt.T[1] = g.TTe; npt.T[2] = g.TTe; @@ -293,7 +293,7 @@ void initializeParticles(SetupParticles& setup_particles, break; case MY_ELECTRON_BG: // electron bg npt.n = g.nb_n0 / sqr(cosh(crd[2] / g.L)); - npt.p[0] = -2. * g.Teb_Te * g.TTe / g.b0 / g.L; + npt.p[0] = 0.;//-2. * g.Teb_Te * g.TTe / g.b0 / g.L; npt.T[0] = g.Teb_Te * g.TTe; npt.T[1] = g.Teb_Te * g.TTe; npt.T[2] = g.Teb_Te * g.TTe;