aethermodel / aether Goto Github PK

Is your feature request related to a problem? Please describe.
The indices system is set up, but need to add a bunch more indices that are needed to run different sub-models

Describe the solution you'd like
Should add:

1. IMF and solar wind (SWMF-type of files & CDAW web files)
2. Ap for MSIS
3. Hemispheric Power
4. AE

Need to add python or other scripts to make these files or commands to download them. The readers should be as close to native format as possible, to minimize preprocessing.

Describe alternatives you've considered
Don't need all of the readers at once. Can implement slowly as needed.

Additional context
A bunch of these readers exist in fortran in GITM. Can rework to C, or can figure out how to get them to work. In the SWMF, this system exists in a util directory, since these are general tools. Maybe we should move this project to a sub-project or something?

Describe the bug
Armadillo was included in include/ folder when that should be strictly Aether stuff. Armadillo is already in share/include/

Expected behavior
No armadillo in include/ folder

Will fix in my next pull request. CMake should recognize armadillo in new folder. If not then that can be a reported bug.

Is your feature request related to a problem? Please describe.

With the addition of the json system, it is clear that users can really mess up the inputs and mistype things, so we need to have some sort of a checker/reporter system for inputs.

Describe the solution you'd like

There should be some sort of reference that has a complete description of variables allowed in the code (e.g., default.json). Then, the modified aether.json file should be checked against this, so that if there are things that are not included that should be, the code can stop, or if there are things that are included that are not described, a warning could be thrown.

Describe alternatives you've considered

The problem with this is that it could be too rigid. Allowing users to mess up allows the code to be extremely flexible, but harder to figure out what could be going wrong with a run that doesn't work. So, I would suggest that the user can throw a flag to check inputs or not.

Additional context

This is a relatively simply thing to implement, so it may be a good first time coding task for an undergrad.

Is your feature request related to a problem? Please describe.

We can use automatic documentation with
Doxygen which is common for C++
specifically.

Describe the solution you'd like

Start documenting code in the .h files with Doxygen

Describe alternatives you've considered

Sphinx, which is popular with Python but I have not seen it used in C++ before.

Is your feature request related to a problem? Please describe.
Aaron's implementation of output.cpp outputs simple netCDF files, while we need to make them a bit more sophisticated so standard netCDF readers and visualizers for atmospheric output files.

Describe the solution you'd like
Look at cfconventions.org for some sort of description.

Describe alternatives you've considered
We could output some other type of data file besides netCDF, but I am not sure what kind.

Additional context
N/A

Is your feature request related to a problem? Please describe.
Need to implement boundary conditions on different variables.

Describe the solution you'd like
Should be in the neutrals and ions classes. Should have (0) reflective (zero at the boundary); (1) fixed value; (2) continuous; (3) continuous gradient; and (4) higher order continuous gradient (?). Type of gradient should be set in the earth.in (planet.in) file for each species and state, maybe? Or in the aether.in file.

Describe alternatives you've considered
Implement each type of BC as a generic system, then call with each variable to get the values. Implement something simple to begin and then make more complex later. All of the coefficients and spatial coefficients should be calculated before hand.

Additional context
GITM has the higher order scheme implemented (thanks Jared!), so we can take this.

Is your feature request related to a problem? Please describe.
In src/bfield.cpp, the following lines state:

else if (input.get_bfield_type() == "dipole") {
    // This is an approximation right now, due to the fact that the
    // pole on some planets (including Earth), have an offset, so the
    // pole location should be altitude dependent.  For many planets,
    // this is a very small error, but should be fixed at some point.

Describe the solution you'd like
I would like this error to be fixed at some point.

Describe alternatives you've considered
Leaving only a simple dipole field

Is your feature request related to a problem? Please describe.
I'm always frustrated when lat and lon are in radians instead of degrees.

Describe the solution you'd like
I would like lat and lon to be outputted in degrees.

Describe alternatives you've considered
An alternative is to leave things in radians and have the conversion to degrees happen in external plotting scripts.

Is your feature request related to a problem? Please describe.
Need to implement viscosity in the momentum equation.

Describe the solution you'd like
Create a file called neutral_momentum_viscosity.cpp and implement in there. In GITM there is the viscosity of the horizontal winds in the vertical direction, but we should implement this on the vertical winds in the horizontal direction. I am not sure, but we should be able to use the conduction solver for the viscosity in the vertical direction.

Describe alternatives you've considered
Should look at GITM and see how it is implemented there. There is no horizontal viscosity, but should be able to implement this with "simple" gradients.

Additional context
Look at GITM for the vertical. Need to implement gradients in the horizontal direction.

Describe the bug

Using a variable such as iErr is not considered best practices. C++ has a
keyword called throw that can send errors that can be descriptive

To Reproduce

throw "Error: File not found"

Expected behavior

catch (const string err_msg) { cerr << err_msg << endl; }

Is your feature request related to a problem? Please describe.
We need parallelization in both ensembles and grid discretization.

Describe the solution you'd like
At first, the parallelization can be without message passing. Need to have each processor read an input file that shows how many processors are devoted to ensembles and blocks (grid). Then, number of lat/lon/alt(?) blocks needs to be read and the code has to figure out grid discretization. Not sure what to do about multiple grids at this point. But, at first, we ignore?

Describe alternatives you've considered
Should have a nompi version. So, can compile with nompi settings. Need preprocessor for this?

Additional context
Nothing.

Is your feature request related to a problem? Please describe.
In order to pass variables to DART, we need to fully describe them.

Describe the solution you'd like
Need to create a specification document and maybe sample code to provide an example. This actually should be done anyways, so that developers can know how to do things.

Describe alternatives you've considered
I am not sure if we should just create some well documented source code to do this, or a document that outlines it. Maybe both.

Additional context
N/A

Is your feature request related to a problem? Please describe.
In the lower portion of the atmosphere, the neutrals collide with each other and their velocities are often tied together, therefore a frictional term needs to exist that binds them together.

Describe the solution you'd like
Take the neutral friction from GITM and implement in Aether.

Describe alternatives you've considered
Implement in the momentum solver in some way? So, solve the coupled equations? I am not sure how this would work. I wouldn't do this immediately.

Additional context
Need the matrix solvers ludcmp and lubksb.

Describe the bug
The netCDF libraries don't seem to work with all compiler / OS combinations, so the code won't compile.

To Reproduce
Need to be on some architecture. For example, it doesn't seem to work with OSX and g++10/11? maybe?

Expected behavior
Code won't compile.

Screenshots
Not going to add these.

System information (please complete the following information):
See above.

Solution

• Have implemented a Config.py script in the main directory to allow configuring the code now.
• Have made two types of outputs - netCDF (output_netcdf.cpp) and binary (output_binary.cpp). Config.py can choose which one is compiled (linked to output.cpp).
• output_binary.cpp outputs a binary file and a header (ascii) file.

Still Needed

• A python reader for the header/binary combo. Should be able to parse data file into variable names, variable units, time, etc, and read in the binary file, creating the proper 3D arrays. This could then be used for plotting if people want.
• A python writer for netCDF. This should take the parsed data that is read in and write out a netCDF that can be read by the normal aetherpy code.
• aetherpy should be modified so that if it gets a ".header" file or a ".bin" file, it looks for the other and reads the appropriate file combo. Remember, GITM uses *.bin files, so it needs to keep working with GITM as well as aether.

Is your feature request related to a problem? Please describe.
Code needs to be adaptive to the current conditions, so dt should be dynamic.

Describe the solution you'd like
Should calculate a dt for neutrals (vertical and horizontal) and ions (vertical and horizontal), then compare and take lowest. Needs to sync across all grid processors but not ensemble processors. (Maybe - may want output for ensembles at the EXACT same time? Maybe implement in the future?)

Describe alternatives you've considered
Straightforward to implement, so shouldn't be a problem.

Additional context
Nope.

Remove doxygen created files and allow the user to make it themselves.
The doxygen files will be in an aether.github.io/docs page

I think static Makefiles can be hard to use. We can keep the static makefile but
I believe we should also have a CMakeLists.txt to use CMake that can handle
dependencies well.

Describe the solution you'd like

Add a CMakeLists.txt

Describe alternatives you've considered

Just using a static Makefile but this might be hard to track all edge cases. And
might be reinventing the wheel.

Additional context

CMake is a well known and well used software. I've used it before and can add it
and maintain it myself.

Is your feature request related to a problem? Please describe.
When modifying output_binary.cpp, you really need to modify output_netcdf.cpp also. This is a pain.

Describe the solution you'd like
Develop a function where you load the arrays (with variable names and such) and then call the output function to dump those variables in either binary or netcdf, disconnecting the different steps. (i.e., call output, which loads the variables into some sort of output structure, and then call dump_output or something that just writes those variables to netcdf or binary.)

Describe alternatives you've considered
Only support one output type, which I don't really like.

Additional context
This would make it incredibly easy for users to make additional output file types.

Is your feature request related to a problem? Please describe.
Some of the code features and attributes could use more explanation

Describe the solution you'd like
Although this is general and ongoing, I would specifically (for this issue) like to see:

• descriptions of the existing classes
• reasons behind different dimensions (nGeoAlts)
• makefile flags

Describe alternatives you've considered
Not making a note of specific areas that I think need comments and leaving it to chance.

Additional context
Items noted for future improvement in PR #31

Is your feature request related to a problem? Please describe.
The science values we care about are all floats float instead of double precision floats double.

Describe the solution you'd like
Change the science values to default to double.

Describe alternatives you've considered
Keeping as float but better precision is better results.

Is your feature request related to a problem? Please describe.
Need to have cooling due to NO, CO2, and O2. These are simple equations (for Earth).

Describe the solution you'd like
Need a planet file that contains planet-specific source terms. This should go into that. Maybe have a separate calculation for each of the different processes (like calc_radiative_NO or something) which can be reused for different planets.

Describe alternatives you've considered
Hack a solution to begin with for Earth, and generalize later when we try to run different planets.

Additional context
None.

Is your feature request related to a problem? Please describe.
The #item type of input that GITM and the SWMF use takes a bit to get used to, while the json format is much more universal. So, it would be good to add json as an input format.

Describe the solution you'd like
Incorporate the json library from (e.g.) https://github.com/nlohmann/json into Aether.

Describe alternatives you've considered
Stick with the #item format for the inputs.

Additional context
Not sure.

Is your feature request related to a problem? Please describe.
We would like the CubeSphere to be out main grid structure for the code. This is for the grid structure in general and not the solvers - that will be a separate issue.

Describe the solution you'd like
There are a bunch of places where we can get a cubesphere grid, but will have to extract it:

• www.gfdl.noaa.gov/fv3
• NRL neptune model (this may be the best place, since it is a non-hydro model, so may be closest to Aether)

Describe alternatives you've considered

• The SWMF has a CubeSphere grid.
• Other generic fluid code solvers also have CubeSphere grids.

Additional context
N/A

Is your feature request related to a problem? Please describe.
Given the electric field, magnetic field, gravity, gradients in pressure, etc, calculate the ion drift (along the field-lines and across, if applicable).

Describe the solution you'd like
This should be a relatively simple code to do these things. Can implement an implicit scheme or a steady-state scheme for the ion velocity. This can be done before the e-field calculation, if needed. The field-aligned component doesn't need this. The question is, which grid will this be done on. For now, do it on the geo grid, I think.

Describe alternatives you've considered
Should be straightforward.

Additional context
Nope.

Is your feature request related to a problem? Please describe.
In the src/electrodynamics_unit_tests.cpp program, there is a section of code that is repeated frequently:

//results verification
  for (int i = 0; i < vals.n_rows; ++i) {
    for (int j = 0; j < vals.n_cols; ++j) {
      if (res(i, j) != answer(i, j))
        return false;
    }
  }

Describe the solution you'd like
This code snippet should be included as an evaluation routine that is called by the test functions.

Describe alternatives you've considered
Leave in duplicate code that makes the routines harder to maintain.

Is your feature request related to a problem? Please describe.
There is a need to have a bunch of rotations and coordinate transformations. Should implement a variety of these and then have unit tests for them.

Describe the solution you'd like
Add to the transform.cpp code to include more rotations (scalar and vector) and coordinate transformations (like GEO->GSE) or (GEO->MAG), etc.

Describe alternatives you've considered
We may not have to implement all of these, since they probably exist all over the web.

Additional context
Nope.

Is your feature request related to a problem? Please describe.
Collisions are critically important when calculating momentum and energy terms, and collision frequencies are important in this. So, we need all of the ion-neutral, ion-ion, ion-electron, and electron-neutral collision terms.

Describe the solution you'd like
A CSV file should have all of the possible collision terms in it. Then the code should read these files and pair them with the existing ions and neutrals.

Describe alternatives you've considered
These are very planet specific, so it is better to have a massive CSV file that can be updated, so we are consistent across the planets. Then, when someone uses a new cross-section, then everyone gets access to this.

Additional context
Need to compare the results to empirical relationships, maybe.

Is your feature request related to a problem? Please describe.
Need chemistry to make the model work.

Describe the solution you'd like
Need to make a generic chemical solver to begin with, where the code takes the sources and losses and calculates the change in density. Sources and losses can be planet specific at this point, but need to migrate to a CSV file with different reaction rates. Ionization rates from the EUV (and eventually aurora). Include heating from chemistry.

Describe alternatives you've considered
At first, set up generic solver for chemistry (sources and losses in each neutral and ion). Implement a hacked simple chemistry scheme for specific reactions, then implement CSV reader for chemistry. The last part will be harder.

Additional context
Pretty critical stuff here.

Is your feature request related to a problem? Please describe.
Vertical velocities of the neutrals will all be different, so we will need to implement a neutral-neutral friction term for the vertical solver.

Describe the solution you'd like
Create a file called neutral_momentum_diffusion.cpp and implement in there. In GITM there is a solver for this (calc_neutral_friction.f90). We should be able to optimize how this is solved given the armadillo library.

Describe alternatives you've considered
This is similar to how the ion - ion collisions should be treated. Right now those are ignored, so we could try to look at solutions that are generalizable for both. Or could hack together directly what is from GITM.

Additional context
This needs to be very "close" to the vertical solver for the neutrals, since the gravity + gradient pressure forces are so large.

Is your feature request related to a problem? Please describe.
After the potential is derived, the electric field needs to be calculated. This is needed to calculate the ion drifts.

Describe the solution you'd like
Take the gradient of the electric potential in all directions. Make sure that the gradient in the field-aligned direction is zero.

Describe alternatives you've considered
Need to make sure that if there is no magnetic field, then the electric field is zero!

Additional context
Nope.

Is your feature request related to a problem? Please describe.
As discussed in the indices->develop pull request, there are a number of variables that are not really named consistently. So, there was a suggestion that we name constants with a bit more consistency.

Describe the solution you'd like
Constants will be named with names that are in CAPS. Constant indices will start with "i" and end with an "". For example, iF107A is a constant index, while PI is a physical constant.

Describe alternatives you've considered
See the linked conversation for more considerations.

Additional context
Nope.

Is your feature request related to a problem? Please describe.
Plotting software would like cell corners. This is sort of important when the grid is not "normal".

Describe the solution you'd like
Should output a single file at the start of the run with complete grid information (cell corners, edges, centers). Then a python code should be written to read and use this grid information. Data should be stored as blocks in individual block files, then a post processing code can make a netcdf of the blocks.

Describe alternatives you've considered
parallel write into netcdf files. Needs to have the right libraries installed.

Additional context
None.

In this example:

missing = 9999.9 was repeated a lot. This can be refactored. There are other parts like this too but I can look around.

Is your feature request related to a problem? Please describe.
The auroral precipitation needs an energy deposition code that takes the spectra and turns it into an altitude profile of heating and ionization rates.

Describe the solution you'd like
First solution could just take the Frahm et al. code and use that (what is in GITM). It is in fortran, so we would have to figure out how to link it.

Describe alternatives you've considered
We should try to get the Fang et al., electron and ion energy deposition code.

Additional context
None.

Is your feature request related to a problem? Please describe.
Ensemble members need to be able to perturb indices or drivers.

Describe the solution you'd like
Need to have a specification of:

• Which index to perturb and how to perturb
• Percent biased perturbations (multiply by % consistently)
• Constant value biased perturbations (+/- value consistently)
• Noise (add ever-changing noise to the system of given standard deviation)

Describe alternatives you've considered
Need a perturbation system, where you feed in a value (or vector of values) and get out perturbed values.

Additional context
N/A

Is your feature request related to a problem? Please describe.
It would be nice if the metadata contained the min and max values for any given output variable.

Describe the solution you'd like
I would like the min and max values of the output data to be calculated and included in the respective metadata.

Describe alternatives you've considered
An alternative is to not include the min and max value.

Additional context
Code seed from output_netcdf.cpp, ln 181:

int min = std::min_element(neutrals.species[iSpecies].density_scgc.begin(), neutrals.species[iSpecies].density_scgc.end());
denVar[iSpecies].putAtt(MIN_VAL, std::min_element(neutrals.species[iSpecies].density_scgc.begin(), neutrals.species[iSpecies].density_scgc.end()));

Is your feature request related to a problem? Please describe.
Joule heating, frictional heating, and heat transfer are key processes and are symmetric in the different equations. Additionally, the momentum transfer terms are symmetric. So, there should be a general place to do these calculations that are agnostic as to the planet.

Describe the solution you'd like
Functions that calculate the general form of the collisional heating and momentum transfer terms. Then specific calls with the ions, neutrals, and electrons.

Describe alternatives you've considered
The real question is whether we will have one ion velocity or ion velocities for each species. Same for temperatures. I don't know how stiff the equations are for the heat transfer, so I don't know whether this is easy to implement or not. For now, I would consider individual velocities and one temperature for bulk ions.

Additional context
Nope.

Is your feature request related to a problem? Please describe.
Need a clear descriptor of how to download, configure, compile, and run the main and "develop" branches of the code.

Describe the solution you'd like
I think that there should be a quickstart in the main directory, maybe? This could point to a more in-depth html version in the documentation.

Describe alternatives you've considered
Could have a PDF file or something that describes stuff - of an html file. I think that we would prefer a text document.

Additional context
N/A

Is your feature request related to a problem? Please describe.
The electron temperature needs to be calculated.

Describe the solution you'd like
Once the collisional heating terms are calculated, there is still a LONG way to go on the electron temperature. This has a HUGE number of terms and is very complicated. We need to calculate this much more carefully than in the past, since the code in GITM is just spaghetti.

Describe alternatives you've considered
Can put in an approximation to begin with.

Additional context
Zhu and Ridley describe the method for the electron temperature equation. GITM's code works.

Is your feature request related to a problem? Please describe.
Some code should exist to get the auroral specification at the top of the model. And there should be some code for calculating the auroral energy deposition. There needs to be code to link the two

Describe the solution you'd like

1. Take auroral patterns and break them into N mono-energetic bins.
2. Integrate the total number density between the current grid and the top of the model. Tilt the grids by the inclination of the field line!
3. Run the energy deposition code to deposit the energy for the different types of aurora.
4. Turn the results into ionization rates and heating rates.

Describe alternatives you've considered
Still need this code to link the two together.

Additional context
Nope.

Is your feature request related to a problem? Please describe.
Implement ion drag in neutral momentum equation - can use collision frequencies from ion velocity calculations

Describe the solution you'd like
Create file called neutrals_momentum_ion_drag.cpp and implement ion drag term.

Describe alternatives you've considered
Above is good.

Additional context
Implemented in GITM, so the forcing term can be copied from there.

We need to make sure to finish the community checklist. I think having a code of
conduct would be good.

We can use a premade code of conduct that github offers.

We may write our own code of conduct. In order to be inclusive I think renaming
our branch main instead of master would be nice and also avoid words like
blacklist and use inclusive alternatives like blocklists

Describe the bug
Aaron created a neutral class with all of the neutrals stored in a C++ vector named "neutral". So, to use something like O, you have use neutral.neutrals[0].density_scgc. The ions name the different constituents "species", so to use something like O+, you have to use ions.species[0].density_scgc. I like "species" better than "neutrals", so all of the "neutrals" should be changed to "species".

To Reproduce
Just look into neutrals.h and ions.h.

Expected behavior
Nothing bad, just stylistically bad.

Screenshots
If applicable, add screenshots to help explain your problem.

System information (please complete the following information):
All

Additional context
Should be easy to fix this. Just will take a bit of time.

Is your feature request related to a problem? Please describe.
Too many 3d loops in the code. It is irritating and there is a need to differentiate the 3d-to-1d look-up between the magnetic and geographic grids, which adds if statements.

Describe the solution you'd like
Need to eliminate loops, and armadillo code does this. Armadillo also takes care of different grids.

Describe alternatives you've considered
Might think of using templates, but I don't know how to do this.

Additional context
Qusai suggested using armadillo, and this is a good solution.

Is your feature request related to a problem? Please describe.
Aether keeps track of time in seconds since Jan 1, 1965. We could store time as a class, and have different methods to get units of time (seconds, UT, array of y, m, d, h, m, s, etc).

Describe the solution you'd like
A class would make this more Class-like, than having just a bunch of functions thrown together to convert to different units.

Describe alternatives you've considered
We could leave this as is, but it would make it much more similar to "datetime" in python to do it this way.

Additional context
N/A

Is your feature request related to a problem? Please describe.
Need to create ensemble members for parallel execution.

Describe the solution you'd like
Need to:

• Create MPI communicator for parallel execution
• Output files that capture ensemble member number
• Output files that capture mean and standard deviation of ensemble (post processing or in code???)

Describe alternatives you've considered
Could have multiple ensemble members on a single processor. I am unsure of how I feel about this. It would be just easier to assume that MPI handles all of this...?

Additional context
N/A

Is your feature request related to a problem? Please describe.
The bulk ion temperature needs to be calculated

Describe the solution you'd like
Once the collisional heating terms are in place, it should be trivial to calculate the ion temperature, since it is mostly a balance between the different collisional heating terms.

Describe alternatives you've considered
The real question is whether to implement a parallel and perpendicular temperature for the ions. Not sure. Need to look at the equations to see how complex they are. This could be put off for a bit.

Additional context
Lindsey whats-her-name from Canada does the parallel and perp temperatures. I need to find her paper.

Describe the bug
The output files are larger and more combersome than they need to be because altitude, latitude, and longitude are included as data variables instead of data coordinates.

To Reproduce
Open a 3DALL file.

Expected behavior
I expect the netcdf file to have coodinates of time, longitude, latitude, and altitude.

System information (please complete the following information):

• OS: OS X
• Compilers used, with version numbers: ncdump/NA
• Other details about your setup that could be relevant

Additional context
This will simplify loading the data into xarray Datasets.

Is your feature request related to a problem? Please describe.
Need a vertical solver for the continuity, momentum, and energy equation.

Describe the solution you'd like
Take the rusanov solver from GITM and implement that first. Need to take all of the state variables and make 1D arrays with them, then solve and remap them back to the 3D solutions.

Describe alternatives you've considered
Once things are more developed, take the AUSM solver and implement that. May want to figure out how we can put in a hydrostatic solver to save time.

Additional context
Nope.

Is your feature request related to a problem? Please describe.
Add a high-latitude potential and aurora to the model. These will be sub-models, either empirical or file driven.

Describe the solution you'd like
In GITM, this is done where you can choose which sub-model you would like to run. You then set the indices and grid and then ask for the potential or aurora. We should use the same framework. If we can get the fortran code to work, then we don't need to do much. If we want to rewrite to C/C++ that is fine. This is not a computational heavy section of the code, so we can clean it up (i.e., update from Fortran into C++) and it will be much nicer.

Describe alternatives you've considered
Figure out how to get the fortran code working and call it good. I already have the sub-model as a stand-alone library, so maybe we just use that to begin with?

Additional context
Some of the codes (e.g., Weimer) are written in Fortran anyways, so we will have to get those working.