@@ -1239,7 +1239,6 @@ class Elm_6f : public PhysicsModel {
12391239 diffusion_coef_Hmode0 = options[" diffusion_coef_Hmode0" ].withDefault (1.0 ); // default value of radial diffusion coefficient
12401240 diffusion_coef_Hmode1 = options[" diffusion_coef_Hmode1" ].withDefault (10.0 ); // upper limit of radial diffusion coefficient
12411241
1242- // path = options["path"].withDefault("./"); // The path of the original Vexb data
12431242 pos_filter = options[" pos_filter" ].withDefault (false ); // switch to turn on the filter of the negative value of zonal background
12441243 pos_filter2 = options[" pos_filter2" ].withDefault (false ); // switch to turn on the filter inside certain position
12451244 pos_filter_zf = options[" pos_filter_zf" ].withDefault (false ); // switch to turn on the filter of the dc profiles inside certain postion with tanh function
@@ -1716,6 +1715,7 @@ class Elm_6f : public PhysicsModel {
17161715 if (mesh->get (Bbar, " bmag" )) { // Typical magnetic field
17171716 Bbar = 1.0 ;
17181717 }
1718+
17191719 // Lbar = options["Lbar"].doc("Reference length scale [m]").withDefault(1.0);
17201720 if (mesh->get (Lbar, " rmag" )) { // Typical length scale
17211721 Lbar = 1.0 ;
@@ -1742,7 +1742,8 @@ class Elm_6f : public PhysicsModel {
17421742
17431743 Tbar = Lbar / Va;
17441744
1745- output.write (" Normalisations: Bbar = {:e} T Lbar = {:e} m\n " , Bbar, Lbar);
1745+ output.write (" Normalisations:\n " );
1746+ output.write (" \t Bbar = {:e} T Lbar = {:e} m\n " , Bbar, Lbar);
17461747 output.write (" \t Va = {:e} m/s Tbar = {:e} s\n " , Va, Tbar);
17471748 output.write (" \t Nbar = {:e} * {:e} m^-3\n " , Nbar, density);
17481749 output.write (" \t Pibar = {:e} Pa Pebar = {:e} Pa\n " , ee * Tibar * Nbar * density, ee * Tebar * Nbar * density);
@@ -1965,6 +1966,7 @@ class Elm_6f : public PhysicsModel {
19651966 P00 = P00 / (ee * Tebar * Nbar * density);
19661967 N0 = n0_height * (pow (P0 / P00, 0.3 ));
19671968 } else {
1969+ // N0 = N0tanh(n0_height, n0_ave, n0_width, n0_center, n0_bottom_x);
19681970 N0 = N0tanh (n0_height * Nbar, n0_ave * Nbar, n0_width, n0_center, n0_bottom_x); // Do we need Nbar here?
19691971 }
19701972 Ti0 = P0 / N0 / (1.0 + Zi);
@@ -2512,11 +2514,12 @@ class Elm_6f : public PhysicsModel {
25122514 // kappa_par_i_lin = B0; // It was B0 in older version of v3
25132515 // kappa_par_e_lin = B0; // It was B0 in older version of v3
25142516
2515- for (jx = 0 ; jx < mesh->LocalNx ; jx++)
2517+ for (jx = 0 ; jx < mesh->LocalNx ; jx++) {
25162518 for (jy = 0 ; jy < mesh->LocalNy ; jy++) {
25172519 kappa_par_i_lin (jx,jy) = kappa_par_i (jx,jy,0 );
25182520 kappa_par_e_lin (jx,jy) = kappa_par_e (jx,jy,0 );
25192521 }
2522+ }
25202523 }
25212524
25222525 if (diffusion_perp > 0.0 ) {
@@ -4563,6 +4566,9 @@ class Elm_6f : public PhysicsModel {
45634566
45644567 if (nonlinear) {
45654568 ddt (U) -= bracket (phi, U, bm_exb); // Advection
4569+ /* if (compress0)
4570+ // ddt(U) -= Vipar*Grad_par(U);
4571+ ddt(U) -= Vpar_Grad_par(Vipar, U);*/
45664572
45674573 if (compress0 && include_vipar) {
45684574 ddt (U) -= Vpar_Grad_par (Vipar, U);
@@ -4808,7 +4814,7 @@ class Elm_6f : public PhysicsModel {
48084814 }
48094815
48104816 if (compress0) {
4811- // ddt(Ni) -= Vipar * Grad_parP(N0, CELL_YLOW );
4817+ // ddt(Ni) -= Vipar * Grad_parP(N0);
48124818 // ddt(Ni) -= Vpar_Grad_par(Vipar, N0);
48134819 if (include_vipar) {
48144820 ddt (Ni) -= Vipar * Grad_parP (N0);
@@ -4835,7 +4841,7 @@ class Elm_6f : public PhysicsModel {
48354841 }
48364842
48374843 if (nonlinear) {
4838- // ddt(Ni) -= Vipar * Grad_par(Ni, CELL_YLOW );
4844+ // ddt(Ni) -= Vipar * Grad_par(Ni);
48394845
48404846 // ddt(Ni) -= Vpar_Grad_par(Vipar, Ni);
48414847 // ddt(Ni) += Vipar * bracket(Psi, N0, bm_mag)*B0;
@@ -4964,7 +4970,7 @@ class Elm_6f : public PhysicsModel {
49644970 }
49654971
49664972 if (nonlinear) {
4967- // ddt(Ti) -= Vipar * Grad_par(Ti, CELL_YLOW );
4973+ // ddt(Ti) -= Vipar * Grad_par(Ti);
49684974
49694975 // ddt(Ti) -= Vpar_Grad_par(Vipar, Ti);
49704976 // ddt(Ti) += Vipar * bracket(Psi, Ti0, bm_mag)*B0;
@@ -4990,48 +4996,47 @@ class Elm_6f : public PhysicsModel {
49904996
49914997 if (diffusion_par > 0.0 ) {
49924998
4993- ddt (Ti) += kappa_par_i * Grad2_par2 (Ti) / N0; // Parallel diffusion
4994- ddt (Ti) += Grad_par (kappa_par_i) * Grad_par (Ti) / N0;
4999+ ddt (Ti) += kappa_par_i * Grad2_par2 (Ti) / N0; // Parallel diffusion
5000+ ddt (Ti) += Grad_par (kappa_par_i) * Grad_par (Ti) / N0;
49955001
4996- if (diff_par_flutter) {
4997- if (nonlinear) {
4998- if (evolve_psi)
4999- bracket1i = -bracket (Psi, Ti + Ti0, bm_mag) * B0;
5000- else
5001- bracket1i = -bracket (Apar, Ti + Ti0, bm_mag);
5002+ if (diff_par_flutter) {
5003+ if (nonlinear) {
5004+ if (evolve_psi)
5005+ bracket1i = -bracket (Psi, Ti + Ti0, bm_mag) * B0;
5006+ else
5007+ bracket1i = -bracket (Apar, Ti + Ti0, bm_mag);
5008+ } else {
5009+ if (evolve_psi)
5010+ bracket1i = -bracket (Psi, Ti0, bm_mag) * B0;
5011+ else
5012+ bracket1i = -bracket (Apar, Ti0, bm_mag);
5013+ }
5014+ mesh->communicate (bracket1i);
5015+ bracket1i.applyBoundary ();
5016+ gradpar_ti = Grad_par (Ti); // NOTE(malamast): Do I need to add Ti0 here?
5017+ mesh->communicate (gradpar_ti);
5018+ gradpar_ti.applyBoundary ();
5019+
5020+ ddt (Ti) += Grad_par (kappa_par_i) * bracket1i / N0;
5021+ ddt (Ti) += kappa_par_i * Grad_par (bracket1i) / N0;
5022+ if (nonlinear) {
5023+ if (evolve_psi) {
5024+ ddt (Ti) -= bracket (Psi, kappa_par_i, bm_mag) * B0 * gradpar_ti / N0;
5025+ ddt (Ti) -= kappa_par_i * bracket (Psi, gradpar_ti, bm_mag) * B0 / N0;
5026+ ddt (Ti) -= bracket (Psi, kappa_par_i, bm_mag) * B0 * bracket1i / N0;
5027+ ddt (Ti) -= kappa_par_i * bracket (Psi, bracket1i, bm_mag) * B0 / N0;
50025028 } else {
5003- if (evolve_psi)
5004- bracket1i = -bracket (Psi, Ti0, bm_mag) * B0;
5005- else
5006- bracket1i = -bracket (Apar, Ti0, bm_mag);
5007- }
5008- mesh->communicate (bracket1i);
5009- bracket1i.applyBoundary ();
5010- gradpar_ti = Grad_par (Ti); // NOTE(malamast): Do I need to add Ti0 here?
5011- mesh->communicate (gradpar_ti);
5012- gradpar_ti.applyBoundary ();
5013-
5014- ddt (Ti) += Grad_par (kappa_par_i) * bracket1i / N0;
5015- ddt (Ti) += kappa_par_i * Grad_par (bracket1i) / N0;
5016- if (nonlinear) {
5017- if (evolve_psi) {
5018- ddt (Ti) -= bracket (Psi, kappa_par_i, bm_mag) * B0 * gradpar_ti / N0;
5019- ddt (Ti) -= kappa_par_i * bracket (Psi, gradpar_ti, bm_mag) * B0 / N0;
5020- ddt (Ti) -= bracket (Psi, kappa_par_i, bm_mag) * B0 * bracket1i / N0;
5021- ddt (Ti) -= kappa_par_i * bracket (Psi, bracket1i, bm_mag) * B0 / N0;
5022- } else {
5023- ddt (Ti) -= bracket (Apar, kappa_par_i, bm_mag) * gradpar_ti / N0;
5024- ddt (Ti) -= kappa_par_i * bracket (Apar, gradpar_ti, bm_mag) / N0;
5025- ddt (Ti) -= bracket (Apar, kappa_par_i, bm_mag) * bracket1i / N0;
5026- ddt (Ti) -= kappa_par_i * bracket (Apar, bracket1i, bm_mag) / N0;
5027- }
5029+ ddt (Ti) -= bracket (Apar, kappa_par_i, bm_mag) * gradpar_ti / N0;
5030+ ddt (Ti) -= kappa_par_i * bracket (Apar, gradpar_ti, bm_mag) / N0;
5031+ ddt (Ti) -= bracket (Apar, kappa_par_i, bm_mag) * bracket1i / N0;
5032+ ddt (Ti) -= kappa_par_i * bracket (Apar, bracket1i, bm_mag) / N0;
50285033 }
5034+ }
50295035
5030- if (output_flux_par) {
5031- heatf_par_flutter_i = -kappa_par_i * bracket1i;
5032- }
5036+ if (output_flux_par) {
5037+ heatf_par_flutter_i = -kappa_par_i * bracket1i;
50335038 }
5034- // }
5039+ }
50355040 }
50365041
50375042 if (diffusion_perp > 0.0 ) {
@@ -5216,49 +5221,48 @@ class Elm_6f : public PhysicsModel {
52165221 }
52175222
52185223 if (diffusion_par > 0.0 ) {
5219- ddt (Te) += kappa_par_e * Grad2_par2 (Te) / Ne0; // Parallel diffusion
5220- ddt (Te) += Grad_par (kappa_par_e) * Grad_par (Te) / Ne0;
5221-
5222- if (diff_par_flutter) {
5223- if (nonlinear) {
5224- if (evolve_psi)
5225- bracket1e = -bracket (Psi, Te + Te0, bm_mag) * B0;
5226- else
5227- bracket1e = -bracket (Apar, Te + Te0, bm_mag);
5224+ ddt (Te) += kappa_par_e * Grad2_par2 (Te) / Ne0; // Parallel diffusion
5225+ ddt (Te) += Grad_par (kappa_par_e) * Grad_par (Te) / Ne0;
5226+
5227+ if (diff_par_flutter) {
5228+ if (nonlinear) {
5229+ if (evolve_psi)
5230+ bracket1e = -bracket (Psi, Te + Te0, bm_mag) * B0;
5231+ else
5232+ bracket1e = -bracket (Apar, Te + Te0, bm_mag);
5233+ } else {
5234+ if (evolve_psi)
5235+ bracket1e = -bracket (Psi, Te0, bm_mag) * B0;
5236+ else
5237+ bracket1e = -bracket (Apar, Te0, bm_mag);
5238+ }
5239+ mesh->communicate (bracket1e);
5240+ bracket1e.applyBoundary ();
5241+ gradpar_te = Grad_par (Te); // NOTE(malamast): Should we add Te0 here?
5242+ mesh->communicate (gradpar_te);
5243+ gradpar_te.applyBoundary ();
5244+
5245+ ddt (Te) += Grad_par (kappa_par_e) * bracket1e / Ne0;
5246+ ddt (Te) += kappa_par_e * Grad_par (bracket1e) / Ne0;
5247+
5248+ if (nonlinear) {
5249+ if (evolve_psi) {
5250+ ddt (Te) -= bracket (Psi, kappa_par_e, bm_mag) * B0 * gradpar_te / Ne0;
5251+ ddt (Te) -= kappa_par_e * bracket (Psi, gradpar_te, bm_mag) * B0 / Ne0;
5252+ ddt (Te) -= bracket (Psi, kappa_par_e, bm_mag) * B0 * bracket1e / Ne0;
5253+ ddt (Te) -= kappa_par_e * bracket (Psi, bracket1e, bm_mag) * B0 / Ne0;
52285254 } else {
5229- if (evolve_psi)
5230- bracket1e = -bracket (Psi, Te0, bm_mag) * B0;
5231- else
5232- bracket1e = -bracket (Apar, Te0, bm_mag);
5233- }
5234- mesh->communicate (bracket1e);
5235- bracket1e.applyBoundary ();
5236- gradpar_te = Grad_par (Te); // NOTE(malamast): Should we add Te0 here?
5237- mesh->communicate (gradpar_te);
5238- gradpar_te.applyBoundary ();
5239-
5240- ddt (Te) += Grad_par (kappa_par_e) * bracket1e / Ne0;
5241- ddt (Te) += kappa_par_e * Grad_par (bracket1e) / Ne0;
5242-
5243- if (nonlinear) {
5244- if (evolve_psi) {
5245- ddt (Te) -= bracket (Psi, kappa_par_e, bm_mag) * B0 * gradpar_te / Ne0;
5246- ddt (Te) -= kappa_par_e * bracket (Psi, gradpar_te, bm_mag) * B0 / Ne0;
5247- ddt (Te) -= bracket (Psi, kappa_par_e, bm_mag) * B0 * bracket1e / Ne0;
5248- ddt (Te) -= kappa_par_e * bracket (Psi, bracket1e, bm_mag) * B0 / Ne0;
5249- } else {
5250- ddt (Te) -= bracket (Apar, kappa_par_e, bm_mag) * gradpar_te / Ne0;
5251- ddt (Te) -= kappa_par_e * bracket (Apar, gradpar_te, bm_mag) / Ne0;
5252- ddt (Te) -= bracket (Apar, kappa_par_e, bm_mag) * bracket1e / Ne0;
5253- ddt (Te) -= kappa_par_e * bracket (Apar, bracket1e, bm_mag) / Ne0;
5254- }
5255+ ddt (Te) -= bracket (Apar, kappa_par_e, bm_mag) * gradpar_te / Ne0;
5256+ ddt (Te) -= kappa_par_e * bracket (Apar, gradpar_te, bm_mag) / Ne0;
5257+ ddt (Te) -= bracket (Apar, kappa_par_e, bm_mag) * bracket1e / Ne0;
5258+ ddt (Te) -= kappa_par_e * bracket (Apar, bracket1e, bm_mag) / Ne0;
52555259 }
5260+ }
52565261
5257- if (output_flux_par) {
5258- heatf_par_flutter_e = -kappa_par_e * bracket1e;
5259- }
5262+ if (output_flux_par) {
5263+ heatf_par_flutter_e = -kappa_par_e * bracket1e;
52605264 }
5261- // }
5265+ }
52625266 }
52635267
52645268 if (diffusion_perp > 0.0 ) {
0 commit comments