Major Update: Expanding bluejay from 1-D to 2-D#11
Conversation
planetarymike
left a comment
planetarymike
left a comment
There was a problem hiding this comment.
Lots of little tweaks requested, but overall this looks great! Thanks for your work to get us to this point.
| # Example modification: Set non-zero flux boundary conditions for O | ||
| # speciesbclist_horiz[:O] = Dict( | ||
| # "f" => [ | ||
| # fill(1e7, num_layers), # Influx at the back edge (cm⁻² s⁻¹) |
There was a problem hiding this comment.
comment here seems to be incorrect if the first entry is for the first / front column edge
There was a problem hiding this comment.
The indexing convention is:
[1] = back edge (first column)
[2] = front edge (last column)
This matches the variable names back_flux and front_flux in the Core.jl implementation, so the comment is correct as written.
There was a problem hiding this comment.
What is the SZA order of the entries in n_current or atm_dict? If I do n_current[:CO2][1], does this give the dayside column (small SZA) or the nightside column?
Maybe your intuition / preference disagrees with mine, but I was thinking that the first column would be the smallest SZA / most dayward, with increasing inxed moving to higher SZA / further to the night side. In this context, I would call the boundary sunward of the first column the "front edge" and the boundary nightward of the last column the "back edge", which is the opposite of what seems to be adopted here.
Interested in @emcangi 's take on this also
There was a problem hiding this comment.
Good point. I updated the comments to make the convention explicit. n_current[sp][1] follows horiz_column_scenario[1], which is the day/smaller-SZA column in the default two-column setup, with increasing ihoriz moving nightward.
The horizontal BC entries are indexed by domain edge: [1] is the lower-index edge outside column 1, and [2] is the higher-index edge outside column n_horiz. I changed the wording to use lower-index/sunward and higher-index/nightward rather than the “front/back”.
There was a problem hiding this comment.
Update: the speciesbclist_horiz edge-flux dictionary (and this example block) was removed entirely along with dx. Horizontal edges are now simply closed (zero flux, the new default) or cyclic, so there is no back/front edge naming left to disambiguate.
The column-ordering part of the question still stands and is now documented in INPUT_PARAMETERS.jl: n_current[sp][1] follows horiz_column_scenario[1], which is the dayside/smallest-SZA column, with increasing index moving nightward – matching your intuition. Transport direction is set by the sign of horiz_transport_timescale_hours (positive = day→night, toward increasing column index).
There was a problem hiding this comment.
all looks ok to me, I guess my intuition is also to have dayside toward the left and nightside toward right
| - Higher-order advection: Upwind schemes are robust but only first-order accurate. Using higher-order methods (e.g., flux limiters, TVD schemes) could reduce numerical diffusion, especially when strong gradients occur across columns. | ||
| - Variable or spatially dependent horizontal winds: The current implementation allows a wind profile per column. More sophisticated models often include spatial variability in the advection term (e.g., using winds from a GCM or prescribing shear with altitude). Extending the wind arrays or reading them from external data would improve realism. | ||
| - Mass-conserving diffusion: The code already averages diffusion coefficients, but ensuring strict global mass conservation (especially with cyclic boundaries) might require verifying that flux leaving one column is exactly balanced by the flux entering the next. Checks or constraints in `fluxcoefs_horiz` could reinforce this. | ||
| - Two-way coupling with vertical transport: Horizontal transport is computed separately from vertical flux coefficients. For atmospheres where horizontal transport interacts with vertical mixing (e.g., along isentropes), a more integrated solver might be beneficial. |
There was a problem hiding this comment.
As long as the entire system ratefn and jacobian is set up to contain both horizontal and vertical transport, the two are solved simultaneously by the code. No operator splitting is used. But maybe this is talking about "diagonal" transport? I'm not sure how our approach would fail under such conditions, but it's possible.
There was a problem hiding this comment.
Yes, the system ratefn and jacobian include both horizontal and vertical transport, and they are solved simultaneously with no operator splitting. The Jacobian in chemJmat_local has entries for above, below (vertical), behind, and infront (horizontal) simultaneously.
There was a problem hiding this comment.
OK, then I think we're OK; situations where horizontal and vertical transport interact should be solvable using the Gear scheme. There may still be issues with our model working in altitude space while the transport physics doesn't really care about altitude, but I would not expect this to be such a big issue that we're unable to obtain scientifically useful results.
|
Thanks so much KP! I updated some comments. |
|
@kpchirakkil I think you may be making edits but I don't see the changes--if you push up your commits, you can then refer to them in your comments like this: "Fixed upwind convection transport in 4989532" (you can just paste the commit hash directly and it will auto shorten to 8 chars) |
# Conflicts: # converge_new_file.jl
emcangi
left a comment
emcangi
left a comment
There was a problem hiding this comment.
I didn't review the plotting routines super closely, but it looks like in all files there's still a lot of instances where nhoriz is not passed as a global variable and I think we mentioned switching it over (and it has been done in some places). I'd recommend Cmd+F'ing for it in all the files to make sure all are passed as global variables
…scales, move horiz_column_width to MODEL_SETUP
…rename n_cur_densities, add docstring, consolidate water enforcement
emcangi
left a comment
emcangi
left a comment
There was a problem hiding this comment.
looks all good to me, only thing i noticed is there are some places where it is declared
n_horiz = GV.n_horiz
which seems unnecessary as GV.n_horiz can just be accessed directly. once Mike approves we can merge I think
| # Horizontal transport timescale in hours. This determines the wind speed via: wind_speed = horiz_column_width / (timescale * 3600) | ||
| # For Venus: 23-44 hours corresponds to 230-120 m/s wind speeds | ||
| # For Mars: Set to 0 for no horizontal transport | ||
| const horiz_transport_timescale = planet == "Venus" ? 30.0 : 0.0 # Baseline shared value (hours) |
There was a problem hiding this comment.
Is anyone else worried about this timescale being in hours when the default time unit throughout the rest of the code is seconds? Am I being too paranoid, or should the unit be carried in the variable name until it's converted?
There was a problem hiding this comment.
Good point – the unit is now carried in the variable name: horiz_transport_timescale_hours (and _hours_neutral / _hours_ion). It's converted to a rate in 1/s in MODEL_SETUP.jl, and only the rate is used downstream.
| # ======================================================================================================= | ||
| const zmin = Dict("Venus"=>90e5, "Mars"=>0.)[planet] | ||
| const dz = 2e5 # Discretized layer thickness | ||
| const horiz_column_width = π * R_P # Surface half-circumference πR (cm) |
There was a problem hiding this comment.
Definitely a quibble, but if the two columns represent average day and average night, then perhaps the representative area-weighted SZA should be 60° for the dayside column and 120° for the nightside, so that the representative horizontal distance is π/3 * R_p .
This doesn't really matter because scaling the distance down would require us to scale the transport timescale down too, in order to get the same model output. Probably this thought can be added as a comment.
There was a problem hiding this comment.
Added the comment in MODEL_SETUP.jl. It keeps horiz_column_width = π * R_P unchanged, but documents the π/3 * R_P point and notes that the chosen width is coupled to the transport timescale used to derive wind speed.
There was a problem hiding this comment.
Follow-up: this is now resolved more thoroughly than by a comment – horiz_column_width (and dx generally) no longer exists in the code, so the comment I mentioned above went away with it. Since the width and the wind speed only ever entered the transport as v/dx = 1/timescale, the distance was redundant with the timescale – scaling one just rescales the other with no effect on the output. The latest push makes that explicit by using the signed timescale as the only input and never forming a velocity or distance. So the π/3 vs π question no longer arises: the timescale is now interpreted directly as the day→night exchange time, whatever representative distance one has in mind.
|
I added some new questions and resolved most of my original comments. I think we are very close to merging, I just want to double-check some of the conventions in conversation before signing off. @emcangi , it looks like there are still a lot of comments from you that could be resolved; I found many of them by expanding the "show more" dialog boxes. |
|
Ok I think I resolved all old stuff tagged to me. But in the process, I found old conversations that @planetarymike started that are unresolved, so now it's Mike's turn to review again! |
3 participants
Expanding bluejay photochemical model from 1-D to 2-D with horizontal coupling. Merged 2-D multicolumn to 1-D master resolving conflicts. Validated the 2-D model for Venus and Mars against 1-D results.