Ionosphere commits #56
Conversation
… and going back into the ionosphere. Warning printed if it happens, could be removed in due time if this doesn't trigger.
…faults to 0 to deactivate and solve every step.
…and a few more [] -> .at().
…at are no longer coupled.
|
Also, I really suggest to do them one by one as e354eac is such a monster. |
sysboundary/ionosphere.cpp
Outdated
|
|
||
| // If we somehow still map into the ionosphere, we missed the 88 degree criterion but shouldn't couple there. | ||
| if(sqrt(x.at(0)*x.at(0) + x.at(1)*x.at(1) + x.at(2)*x.at(2)) < Ionosphere::innerRadius) { | ||
| # TODO drop this warning if it never occurs? To be followed. |
There was a problem hiding this comment.
You, sir, are writing too much python. This is not a c++ source code comment. :)
There was a problem hiding this comment.
Hupsis. I wanted to make that a #warning, or was it #pragma warning? I'll fix it.
|
The last two are waiting to run on Mahti. |
…e-commits Conflicts: sysboundary/ionosphere.cpp
sysboundary/ionosphere.cpp
Outdated
| // If we somehow still map into the ionosphere, we missed the 88 degree criterion but shouldn't couple there. | ||
| if(sqrt(x.at(0)*x.at(0) + x.at(1)*x.at(1) + x.at(2)*x.at(2)) < Ionosphere::innerRadius) { | ||
| # TODO drop this warning if it never occurs? To be followed. | ||
| cerr << "Triggered mapping back into Earth\n"; |
There was a problem hiding this comment.
To my surprise, I see this warning being printed quite often in a small testrun I just set up in my fridge.
There was a problem hiding this comment.
In a small test, with coarse resolution, I can see how it would just miss the threshold and step past the 88° limit.
There was a problem hiding this comment.
We can make it more complex and look e.g. at passing an extremum/return point but the more complex we make it, the more I fear we'll have problems once we have non-trivial B configurations.
…e-commits Conflicts: sysboundary/ionosphere.cpp
…e-commits Conflicts: sysboundary/ionosphere.cpp
…e-commits Conflicts: datareduction/datareducer.cpp
|
I added one typo fix but more usefully that more detailed reporting to logfile instead of the slurm file. Example output:
I thought I wouldn't explicitly add converged/not converged but can be added if requested. |
|
So the last two I guess are uncontroversial (the latter was tedious and is useful for debugging, so unless you have a more compact way of writing that is understandable, it would be nice to keep it I think). |
Density and temperature are the average of nearest neighbors while velocity comes from the ionospheric potential.
Include tstep and time on line for easier grepping and plotting. Separate potential min max difference for N and S hemispheres.
|
Comments on bdbbff0? This makes life easier when grepping/plotting from the logfile. |
|
Of course I forgot f0450bb that gives the separate N and S potential min and max. |
|
21f43b6 fixes the coupling interval feature so that the Vlasov boundary VDFs are also only updated at the set interval like the ionosphere. The if condition in the vlasovBoundary function is a bit ugly due to the counter being incremented in main(), feel free to suggest a better solution if this hurts your eye too much. |
|
5ad350f is a bit WIP. Let me know if this is good enough for now, then I'll clean up the commented regions. If you have another idea, please post. Of course we can also ponder whether we want to find which derivatives to zero out and which to keep, but I'm not sure we want to go there... |
|
Alrighty, thumbs up to the whole bunch of recent commits there! But I cherry-picked them over nonetheless. :) |
…e-commits Conflicts: sysboundary/ionosphere.cpp vlasiator.cpp
| } | ||
|
|
||
| // Calculate Hall and Pedersen conductivity coefficient based on charge carrier density | ||
| const Real Bval = 5e-5; // TODO: Hardcoded B strength here? | ||
| const Real NO_gyroFreq = physicalconstants::CHARGE * Bval / (31*physicalconstants::MASS_PROTON); // Ion (NO+) gyration frequency |
There was a problem hiding this comment.
Here's one diff not due to my .at() craze.
There was a problem hiding this comment.
However, gumics also uses 31 here. I'm still not sure if that is not just simply wrong.
| Real halfdx = 1000 * 0.5 * (atmosphere[h].altitude - atmosphere[h-1].altitude); | ||
| Real halfCH = halfdx * 0.5 * (atmosphere[h-1].hallcoeff + atmosphere[h].hallcoeff); | ||
| Real halfCP = halfdx * 0.5 * (atmosphere[h-1].pedersencoeff + atmosphere[h].pedersencoeff); | ||
| Real halfCpara = halfdx * 0.5 * (atmosphere[h-1].parallelcoeff + atmosphere[h].parallelcoeff); | ||
| Real halfdx = 1000 * 0.5 * (atmosphere.at(h).altitude - atmosphere.at(h-1).altitude); | ||
|
|
||
| nodes[n].parameters[ionosphereParameters::SIGMAP] += (electronDensity[h]+electronDensity[h-1]) * halfCP; | ||
| nodes[n].parameters[ionosphereParameters::SIGMAH] += (electronDensity[h]+electronDensity[h-1]) * halfCH; | ||
| nodes[n].parameters[ionosphereParameters::SIGMAPARALLEL] += (electronDensity[h]+electronDensity[h-1]) * halfCpara; | ||
| nodes.at(n).parameters.at(ionosphereParameters::SIGMAP) += halfdx * 0.5 * (electronDensity.at(h)*atmosphere.at(h).pedersencoeff + electronDensity.at(h-1)*atmosphere.at(h-1).pedersencoeff); | ||
| nodes.at(n).parameters.at(ionosphereParameters::SIGMAH) += halfdx * 0.5 * (electronDensity.at(h)*atmosphere.at(h).hallcoeff + electronDensity.at(h-1)*atmosphere.at(h-1).hallcoeff); | ||
| nodes.at(n).parameters.at(ionosphereParameters::SIGMAPARALLEL) += halfdx * 0.5 * (electronDensity.at(h)*atmosphere.at(h).parallelcoeff + electronDensity.at(h-1)*atmosphere.at(h-1).parallelcoeff); |
There was a problem hiding this comment.
This is also still open.
Conflicts: ioread.cpp
Here's a string of commits I made of all my debugging and developments of yesterday. I suggest cherry-picking one by one.
e354eac was so much work I'd like to have it in. Makes it safer, really...
5462ca0 is cosmetic, but doesn't harm either.
a62574f is the big one.
3333680 is probably better to have in.
56e0020 and 39498b2 are features from gumics, might be useful for further testing and benchmarking.