Skip to content

Commit

Permalink
Merge pull request #145 from AetherModel/IEaddition
Browse files Browse the repository at this point in the history 
Add Fortran IE models such as Weimer and FTA
  • Loading branch information
@aaronjridley
aaronjridley authored Nov 21, 2023
2 parents 3f931fe + 0003b19 commit e986c88
Show file tree
Hide file tree
Showing 13 changed files with 219 additions and 111 deletions.
3 changes: 0 additions & 3 deletions include/neutrals.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -178,9 +178,6 @@ class Neutrals {
/// Vector of all species-specific items:
std::vector<species_chars> species;

/// when computing dt, derive a dt for neutrals:
precision_t dt;

/// Maximum Chapman integral (will give nearly infinite tau in EUV)
precision_t max_chapman = 1.0e26;

Expand Down
7 changes: 7 additions & 0 deletions include/tools.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,13 @@ struct mat_2x2{
arma_mat A22;
};

// ----------------------------------------------------------------------------
// Fix corners in an arma cube
// - basically fill in the corners with values near them
// ----------------------------------------------------------------------------

void fill_corners(arma_cube &values, int64_t nGCs);

// -----------------------------------------------------------------------------
// add cMember into a string just before last period
// -----------------------------------------------------------------------------
Expand Down
17 changes: 13 additions & 4 deletions src/advance.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,9 @@ bool advance(Planets &planet,
report.print(-1, "(1) What function is this " +
input.get_student_name() + "?");

if (didWork & input.get_check_for_nans())
didWork = neutrals.check_for_nonfinites();

gGrid.calc_sza(planet, time);
neutrals.calc_mass_density();
neutrals.calc_mean_major_mass();
Expand All @@ -56,11 +59,17 @@ bool advance(Planets &planet,
// first

neutrals.calc_scale_height(gGrid);
didWork = neutrals.set_bcs(gGrid, time, indices);

if (didWork)
didWork = neutrals.set_bcs(gGrid, time, indices);

if (input.get_nAltsGeo() > 1)
neutrals.advect_vertical(gGrid, time);

if (didWork & input.get_check_for_nans())
didWork = neutrals.check_for_nonfinites();


// ------------------------------------
// Calculate source terms next:

Expand Down Expand Up @@ -97,9 +106,6 @@ bool advance(Planets &planet,
if (input.get_NO_cooling())
neutrals.calc_NO_cool();

neutrals.calc_conduction(gGrid, time);
chemistry.calc_chemistry(neutrals, ions, time, gGrid);

neutrals.vertical_momentum_eddy(gGrid);
calc_ion_collisions(neutrals, ions);
calc_neutral_friction(neutrals);
Expand All @@ -124,6 +130,9 @@ bool advance(Planets &planet,
}
}

if (didWork & input.get_check_for_nans())
didWork = neutrals.check_for_nonfinites();

if (didWork)
didWork = output(neutrals, ions, gGrid, time, planet);

Expand Down
2 changes: 1 addition & 1 deletion src/calc_chemical_sources.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,7 @@ void Chemistry::calc_chemical_sources(Neutrals &neutrals,

for (iReaction = 0; iReaction < nReactions; iReaction++) {

if (report.test_verbose(3)) {
if (report.test_verbose(4)) {
std::cout << "===> Reaction : " << iReaction
<< " of " << nReactions << "\n";
display_reaction(reactions[iReaction]);
Expand Down
22 changes: 16 additions & 6 deletions src/calc_neutral_derived.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -337,6 +337,8 @@ precision_t Neutrals::calc_dt(Grid grid) {
static int iFunction = -1;
report.enter(function, iFunction);

precision_t dt;

if (input.get_is_cubesphere())
dt = calc_dt_cubesphere(grid);
else {
Expand Down Expand Up @@ -382,6 +384,7 @@ precision_t Neutrals::calc_dt_cubesphere(Grid grid) {
static int iFunction = -1;
report.enter(function, iFunction);

precision_t dt;
int iDir;

arma_vec dta(4);
Expand All @@ -390,7 +393,7 @@ precision_t Neutrals::calc_dt_cubesphere(Grid grid) {
int64_t nAlts = grid.get_nAlts();
int64_t nXs = grid.get_nLons();
int64_t nYs = grid.get_nLats();

// dtx dty for reference coordinate system
arma_cube dtx(size(cMax_vcgc[0]));
arma_cube dty(size(cMax_vcgc[0]));
Expand All @@ -401,11 +404,18 @@ precision_t Neutrals::calc_dt_cubesphere(Grid grid) {
// Loop through altitudes
for (int iAlt = 0; iAlt < nAlts; iAlt++) {
// Conver cMax to contravariant velocity first
arma_mat u1 = cMax_vcgc[0].slice(iAlt) % grid.A11_inv_scgc.slice(iAlt) + cMax_vcgc[1].slice(iAlt) % grid.A12_inv_scgc.slice(iAlt);
arma_mat u2 = cMax_vcgc[0].slice(iAlt) % grid.A21_inv_scgc.slice(iAlt) + cMax_vcgc[1].slice(iAlt) % grid.A22_inv_scgc.slice(iAlt);

dtx.slice(iAlt) = grid.drefx(iAlt)*dummy_1 / u1;
dty.slice(iAlt) = grid.drefy(iAlt)*dummy_1 / u2;
arma_mat u1 = sqrt(
cMax_vcgc[0].slice(iAlt) % grid.A11_inv_scgc.slice(iAlt) %
cMax_vcgc[0].slice(iAlt) % grid.A11_inv_scgc.slice(iAlt) +
cMax_vcgc[1].slice(iAlt) % grid.A12_inv_scgc.slice(iAlt) %
cMax_vcgc[1].slice(iAlt) % grid.A12_inv_scgc.slice(iAlt));
arma_mat u2 = sqrt(
cMax_vcgc[0].slice(iAlt) % grid.A21_inv_scgc.slice(iAlt) %
cMax_vcgc[0].slice(iAlt) % grid.A21_inv_scgc.slice(iAlt) +
cMax_vcgc[1].slice(iAlt) % grid.A22_inv_scgc.slice(iAlt) %
cMax_vcgc[1].slice(iAlt) % grid.A22_inv_scgc.slice(iAlt));
dtx.slice(iAlt) = grid.drefx(iAlt) * dummy_1 / u1;
dty.slice(iAlt) = grid.drefy(iAlt) * dummy_1 / u2;
}

// simply some things, and just take the bulk value for now:
Expand Down
151 changes: 79 additions & 72 deletions src/electrodynamics.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@


// -----------------------------------------------------------------------------
// Call (fortran) ionospheric electrodynamics if we enable this
// Call (fortran) ionospheric electrodynamics if we enable this
// in the cmake file
// -----------------------------------------------------------------------------

Expand All @@ -27,36 +27,40 @@ Electrodynamics::Electrodynamics(Times time) {
HaveFortranIe = false;
isSet = false;

#ifdef FORTRAN
// Initialize IE components (reading in the data files):
std::string ieDir = input.get_electrodynamics_dir();
int* ieDir_iArray = copy_string_to_int(ieDir);
std::string efield = input.get_potential_model();
std::string aurora = input.get_diffuse_auroral_model();
#ifdef FORTRAN
// Initialize IE components (reading in the data files):
std::string ieDir = input.get_electrodynamics_dir();
int* ieDir_iArray = copy_string_to_int(ieDir);
std::string efield = input.get_potential_model();
std::string aurora = input.get_diffuse_auroral_model();

if (efield.length() == 0 & aurora.length() == 0) {
HaveFortranIe = false;
if (efield.length() == 0 & aurora.length() == 0)
HaveFortranIe = false;

else {
int* efield_iArray = copy_string_to_int(efield);
int* aurora_iArray = copy_string_to_int(aurora);
ie_init_library(ieDir_iArray, efield_iArray, aurora_iArray, &iError);

if (iError > 0) {
report.print(0, "Error in setting fortran IE!");
IsOk = false;
} else {
int* efield_iArray = copy_string_to_int(efield);
int* aurora_iArray = copy_string_to_int(aurora);
ie_init_library(ieDir_iArray, efield_iArray, aurora_iArray, &iError);
if (iError > 0) {
report.print(0,"Error in setting fortran IE!");
IsOk = false;
} else {
HaveFortranIe = true;
isSet = true;
}
HaveFortranIe = true;
isSet = true;
}
isCompiled = true;
#else
isCompiled = false;
#endif
}

isCompiled = true;
#else
isCompiled = false;
#endif

// If we don't set IE through Fortran, then try to set it through a file:

if (!HaveFortranIe) {
std::string electrodynamics_file = input.get_electrodynamics_file();

if (electrodynamics_file.length() > 0 &
electrodynamics_file != "none") {
// This function sets HaveElectrodynamicsFile = true.
Expand All @@ -83,7 +87,7 @@ Electrodynamics::Electrodynamics(Times time) {
// -----------------------------------------------------------------------------

void Electrodynamics::set_all_indices_for_ie(Times time,
Indices &indices) {
Indices &indices) {
std::string function = "Electrodynamics::set_all_indices_for_ie";
static int iFunction = -1;
report.enter(function, iFunction);
Expand Down Expand Up @@ -111,17 +115,18 @@ void Electrodynamics::set_all_indices_for_ie(Times time,
std::cout << "sw v : " << iVx_ << " " << swv << "\n";
std::cout << "sw n : " << iN_ << " " << swn << "\n";
}
#ifdef FORTRAN
ie_set_time(&time_now);
ie_set_imfby(&imfby);
ie_set_imfbz(&imfbz);
ie_set_swv(&swv);
ie_set_swn(&swn);
ie_set_ae(&ae);
ie_set_au(&au);
ie_set_al(&al);
ie_set_hp_from_ae(&ae);
#endif

#ifdef FORTRAN
ie_set_time(&time_now);
ie_set_imfby(&imfby);
ie_set_imfbz(&imfbz);
ie_set_swv(&swv);
ie_set_swn(&swn);
ie_set_ae(&ae);
ie_set_au(&au);
ie_set_al(&al);
ie_set_hp_from_ae(&ae);
#endif

report.exit(function);
return;
Expand Down Expand Up @@ -154,60 +159,62 @@ bool Electrodynamics::update(Planets planet,
ions.eflux.zeros();
ions.avee.ones();

#ifdef FORTRAN
if (HaveFortranIe) {
report.print(3, "Using Fortran Electrodynamics!");
set_all_indices_for_ie(time, indices);

if (!IsAllocated) {
int nXs = gGrid.get_nX();
ie_set_nxs(&nXs);
int nYs = gGrid.get_nY();
ie_set_nys(&nYs);
int64_t iTotal = nXs * nYs;
mlt2d = static_cast<float*>(malloc(iTotal * sizeof(float)));
lat2d = static_cast<float*>(malloc(iTotal * sizeof(float)));
pot2d = static_cast<float*>(malloc(iTotal * sizeof(float)));
eflux2d = static_cast<float*>(malloc(iTotal * sizeof(float)));
avee2d = static_cast<float*>(malloc(iTotal * sizeof(float)));
IsAllocated = true;
}
#ifdef FORTRAN

if (HaveFortranIe) {
report.print(3, "Using Fortran Electrodynamics!");
set_all_indices_for_ie(time, indices);

if (!IsAllocated) {
int nXs = gGrid.get_nX();
ie_set_nxs(&nXs);
int nYs = gGrid.get_nY();
ie_set_nys(&nYs);
int64_t iTotal = nXs * nYs;
mlt2d = static_cast<float*>(malloc(iTotal * sizeof(float)));
lat2d = static_cast<float*>(malloc(iTotal * sizeof(float)));
pot2d = static_cast<float*>(malloc(iTotal * sizeof(float)));
eflux2d = static_cast<float*>(malloc(iTotal * sizeof(float)));
avee2d = static_cast<float*>(malloc(iTotal * sizeof(float)));
IsAllocated = true;
}

int64_t nZs = gGrid.get_nZ();
int64_t iZ;
int64_t nZs = gGrid.get_nZ();
int64_t iZ;

int iError;
int iError;

for (iZ = 0; iZ < nZs; iZ++) {
copy_mat_to_array(gGrid.magLocalTime_scgc.slice(iZ), mlt2d, true);
copy_mat_to_array(gGrid.magLat_scgc.slice(iZ), lat2d, true);
for (iZ = 0; iZ < nZs; iZ++) {
copy_mat_to_array(gGrid.magLocalTime_scgc.slice(iZ), mlt2d, true);
copy_mat_to_array(gGrid.magLat_scgc.slice(iZ), lat2d, true);

ie_set_mlts(mlt2d, &iError);
ie_set_lats(lat2d, &iError);
ie_update_grid(&iError);
ie_set_mlts(mlt2d, &iError);
ie_set_lats(lat2d, &iError);
ie_update_grid(&iError);

ie_get_potential(pot2d, &iError);
copy_array_to_mat(pot2d, ions.potential_scgc.slice(iZ), true);
ie_get_potential(pot2d, &iError);
copy_array_to_mat(pot2d, ions.potential_scgc.slice(iZ), true);

if (iZ == nZs-1) {
ie_get_electron_diffuse_aurora(eflux2d, avee2d, &iError);
copy_array_to_mat(avee2d, ions.avee, true);
copy_array_to_mat(eflux2d, ions.eflux, true);
}
if (iZ == nZs - 1) {
ie_get_electron_diffuse_aurora(eflux2d, avee2d, &iError);
copy_array_to_mat(avee2d, ions.avee, true);
copy_array_to_mat(eflux2d, ions.eflux, true);
}
}
#endif
}

#endif

if (HaveElectrodynamicsFile) {
report.print(3, "Setting electrodynamics from file!");
auto electrodynamics_values =
get_electrodynamics(gGrid.magLat_scgc,
gGrid.magLocalTime_scgc);
gGrid.magLocalTime_scgc);
ions.potential_scgc = std::get<0>(electrodynamics_values);
ions.eflux = std::get<1>(electrodynamics_values);
ions.avee = std::get<2>(electrodynamics_values);
}
}
}

report.exit(function);
return true;
Expand Down
11 changes: 11 additions & 0 deletions src/exchange_messages_v2.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -601,16 +601,27 @@ void Grid::exchange(arma_cube &data, const bool pole_inverse) {

bool Neutrals::exchange(Grid &grid) {
// For each species, exchange if its DoAdvect is true
int64_t nGCs = grid.get_nGCs();

for (int i = 0; i < nSpecies; ++i) {
if (species[i].DoAdvect)
grid.exchange(species[i].density_scgc, false);

fill_corners(species[i].density_scgc, nGCs);
}

// Exchange temperature
grid.exchange(temperature_scgc, false);
// Fix corners:
fill_corners(temperature_scgc, nGCs);

// Exchange velocity
grid.exchange(velocity_vcgc[0], true);
grid.exchange(velocity_vcgc[1], true);
grid.exchange(velocity_vcgc[2], false);

for (int iDir = 0; iDir < 3; iDir++)
fill_corners(velocity_vcgc[iDir], nGCs);

return true;
}
11 changes: 6 additions & 5 deletions src/init_geo_grid.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -984,6 +984,12 @@ bool Grid::init_geo_grid(Quadtree quadtree,

// Calculate the radius (for spherical or non-spherical)
fill_grid_radius(planet);

// Correct the reference grid with correct length scale:
// (with R = actual radius)
if (input.get_is_cubesphere())
correct_xy_grid(planet);

// Calculate grid spacing
calc_grid_spacing(planet);
//calculate radial unit vector (for spherical or oblate planet)
Expand All @@ -994,11 +1000,6 @@ bool Grid::init_geo_grid(Quadtree quadtree,
// Calculate magnetic field and magnetic coordinates:
fill_grid_bfield(planet);

// Correct the reference grid with correct length scale:
// (with R = actual radius)
if (input.get_is_cubesphere())
correct_xy_grid(planet);

// Throw a little message for students:
report.student_checker_function_name(input.get_is_student(),
input.get_student_name(),
Expand Down
Loading

0 comments on commit e986c88

Please sign in to comment.