Skip to content

Commit

Permalink
FEAT: implement bulk temperature calculation
Browse files Browse the repository at this point in the history
  • Loading branch information
@aaronjridley
aaronjridley committed Nov 22, 2024
1 parent c3d04d0 commit f72443a
Showing 1 changed file with 47 additions and 58 deletions.
105 changes: 47 additions & 58 deletions src/calc_ion_temperature.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -78,79 +78,79 @@ void Ions::calc_ion_temperature(Neutrals neutrals, Grid grid,
// // First ion species only, currently is O+
// nSpecs = 1;
// else
nSpecs = nSpecies;
nSpecs = nSpecies;

if (report.test_verbose(4)) {
std::cout << "Bulk ion temp flag: " << input.get_do_calc_bulk_ion_temp() ?
"true" : "false";
std::cout << " so 'number of ions' is " << nSpecs << "\n";
}

// Loop over all species or assume only bulk calculation
for (iIon = 0; iIon < nSpecs; iIon++) {
//std::cout << "iIon : " << iIon << "\n";
Mi = species[iIon].mass / cAMU;
calc_lambda();

if (input.get_do_calc_bulk_ion_temp()) {
for (iLon = nGCs; iLon < nLons - nGCs; iLon++) {
for (iLat = nGCs; iLat < nLats - nGCs; iLat++) {
//std::cout << "iLon, iLat : " << iLon << " " << iLat << "\n";
// ---------------------------------------------------------------------
// Calculate heat flux (conduction) in 1D; loop over all lat,lon
// ---------------------------------------------------------------------
ratios.zeros();
for (jIon = 0; jIon < nSpecies; jIon++) {
if (jIon != iIon) {
Mj = species[jIon].mass / cAMU;
density_ratio = species[jIon].density_scgc.tube(iLon, iLat) /
species[iIon].density_scgc.tube(iLon, iLat);
density_ratio.clamp(0.001, 1000.0);
ratios = ratios + density_ratio *
(species[jIon].charge * species[jIon].charge /
species[iIon].charge / species[iIon].charge) *
sqrt(Mj / (Mi + Mj)) *
(3 * Mi * Mi + 1.6 * Mi * Mj + 1.3 * Mj * Mj)/
((Mi + Mj) * (Mi + Mj));
}
}
temp1d = species[iIon].temperature_scgc.tube(iLon, iLat);
temp1d(0) = neutrals.temperature_scgc(iLon, iLat, 0);
temp1d(1) = neutrals.temperature_scgc(iLon, iLat, 1);
lambda1d = 3.1e6 / sqrt(Mi) / pow(species[iIon].charge, 4) *
pow(temp1d, 2.5) % (1 + 1.75 * ratios) * cE;
temp1d = temperature_scgc.tube(iLon, iLat);
lambda1d = lambda.tube(iLon, iLat);
lambda1d(1) = lambda1d(2);
lambda1d(0) = lambda1d(2);
//lambda1d = 25.0 * cKB * pow(temp1d, 2.5) * (cKB / species[iIon].mass)
// / species[iIon].nu_ion_ion[iIon] / 8.0;
front1d = 3.0 / 2.0 * cKB * species[iIon].density_scgc.tube(iLon, iLat);
front1d = 3.0 / 2.0 * cKB * density_scgc.tube(iLon, iLat);
dalt1d = grid.dalt_lower_scgc.tube(iLon, iLat);
sources1d = (heating_neutral_friction_scgc.tube(iLon, iLat) +
heating_neutral_heat_transfer_scgc.tube(iLon, iLat));
sources1d = sources1d / front1d;

//std::cout << "iIon lab : " << iIon << "\n" << lambda1d << "\n source:\n " << sources1d << "\n";

conduction1d.zeros(); // reset temp variable to zero

conduction1d = solver_conduction(temp1d,
lambda1d,
front1d,
sources1d,
dalt1d,
dt/100.,
dt/10.,
nGCs,
false);

// The conduction solver gives Tnew-Told, so divide by dt
conduction1d.clamp(200, 5000);
species[iIon].temperature_scgc.tube(iLon, iLat) = conduction1d;

//std::cout << "temp : " << conduction1d << "\n";

} // Lats
} // Lons
} // Ions
temperature_scgc.tube(iLon, iLat) = conduction1d;
}
}
for (iIon = 0; iIon < nSpecies; iIon++)
species[iIon].temperature_scgc = temperature_scgc;

//if (!input.get_do_calc_bulk_ion_temp()) {
// Use the density averaged temperature to fill the bulk temperature
} else {

// Loop over all species or assume only bulk calculation
for (iIon = 0; iIon < nSpecs; iIon++) {
for (iLon = nGCs; iLon < nLons - nGCs; iLon++) {
for (iLat = nGCs; iLat < nLats - nGCs; iLat++) {
temp1d = species[iIon].temperature_scgc.tube(iLon, iLat);
temp1d(0) = neutrals.temperature_scgc(iLon, iLat, 0);
temp1d(1) = neutrals.temperature_scgc(iLon, iLat, 1);
lambda1d = species[iIon].lambda.tube(iLon, iLat);
lambda1d(1) = lambda1d(2);
lambda1d(0) = lambda1d(2);
front1d = 3.0 / 2.0 * cKB * species[iIon].density_scgc.tube(iLon, iLat);
dalt1d = grid.dalt_lower_scgc.tube(iLon, iLat);
sources1d = (species[iIon].heating_neutral_friction_scgc.tube(iLon, iLat) +
species[iIon].heating_neutral_heat_transfer_scgc.tube(iLon, iLat));
sources1d = sources1d / front1d;

conduction1d.zeros(); // reset temp variable to zero
conduction1d = solver_conduction(temp1d,
lambda1d,
front1d,
sources1d,
dalt1d,
dt/10.,
nGCs,
false);

// The conduction solver gives Tnew-Told, so divide by dt
conduction1d.clamp(200, 5000);
species[iIon].temperature_scgc.tube(iLon, iLat) = conduction1d;
} // Lats
} // Lons
} // Ions
tempT.zeros();
tempD.zeros();

Expand All @@ -162,19 +162,8 @@ void Ions::calc_ion_temperature(Neutrals neutrals, Grid grid,
}

temperature_scgc = tempT / tempD;
//}
/*
if (input.get_do_calc_bulk_ion_temp()) {
// Add temperature terms together to advance bulk ion temperature
temperature_scgc = temperature_scgc + dt * (conduction_scgc);
// Use the bulk ion temperature to fill all ion specie temperatures
for (iIon = 0; iIon < nSpecies; iIon++)
species[iIon].temperature_scgc = temperature_scgc;
}
*/

//std::cout << "ion temp : " << temperature_scgc(2,2,5) << " " << neutrals.temperature_scgc(2,2,5) << "\n";
report.exit(function);
return;
}

0 comments on commit f72443a

Please sign in to comment.