Skip to content

Commit 94fcdea

Browse files
abukowski21abukowski21
committed
FEAT: Get vertical neutral (lower) bc's working with new dipole grid
1 parent 6bb4d6d commit 94fcdea

File tree

3 files changed

+98
-85
lines changed

3 files changed

+98
-85
lines changed

‎src/advance.cpp

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -106,8 +106,8 @@ bool advance(Planets &planet,
106106
if (didWork)
107107
didWork = ions.set_bcs(gGrid, time, indices);
108108

109-
//if (didWork)
110-
// didWork = neutralsMag.set_bcs(mGrid, time, indices);
109+
if (didWork)
110+
didWork = neutralsMag.set_bcs(mGrid, time, indices);
111111

112112
didWork = neutralsMag.check_for_nonfinites("Ion Grid: set bcs");
113113

‎src/neutrals_bcs.cpp

Lines changed: 84 additions & 73 deletions
Original file line numberDiff line numberDiff line change
@@ -92,8 +92,8 @@ bool Neutrals::set_upper_bcs(Grid grid) {
9292

9393
h = species[iSpecies].scale_height_scgc.slice(iAlt);
9494
species[iSpecies].density_scgc.slice(iAlt) =
95-
species[iSpecies].density_scgc.slice(iAlt - 1) %
96-
exp(-grid.dk_edge_m.slice(iAlt) / h);
95+
species[iSpecies].density_scgc.slice(iAlt - 1) %
96+
exp(-grid.dk_edge_m.slice(iAlt) / h);
9797
}
9898
}
9999

@@ -118,14 +118,19 @@ bool Neutrals::set_lower_bcs(Grid grid,
118118
json bcs = input.get_boundary_condition_types();
119119
int64_t nGCs = grid.get_nGCs();
120120
int64_t iSpecies, iAlt, iDir;
121+
int64_t nLats = grid.get_nLats();
122+
int64_t nLons = grid.get_nLons();
121123

122124
std::string bcsType = mklower(bcs["type"]);
123125

124126
//-----------------------------------------------
125127
// MSIS BCs - only works if FORTRAN is enabled!
126128
//-----------------------------------------------
127129

128-
if (bcsType == "msis") {
130+
// ALB changes to lower BCs only really work now for dipole grid. Don't use msis
131+
// if we are handed a dipole grid.
132+
133+
if (bcsType == "msis" && !grid.IsDipole) {
129134

130135
report.print(2, "Using MSIS for Boundary Conditions");
131136

@@ -168,78 +173,84 @@ bool Neutrals::set_lower_bcs(Grid grid,
168173
std::cout << " Found in MSIS!\n";
169174

170175
species[iSpecies].density_scgc.slice(0) =
171-
msis.get_mat(species[iSpecies].cName);
176+
msis.get_mat(species[iSpecies].cName);
172177
} else {
173178
if (report.test_verbose(3))
174179
std::cout << " NOT Found in MSIS - setting constant\n";
175180

176181
species[iSpecies].density_scgc.slice(0).
177-
fill(species[iSpecies].lower_bc_density);
182+
fill(species[iSpecies].lower_bc_density);
178183
}
179184

180185
}
181186

182187
} // type == Msis
183188

189+
precision_t sh_ave;
190+
184191
//-----------------------------------------------
185-
// Planet BCs - set to fixed constant values.
192+
// Fill the lower+ ghost cells
186193
//-----------------------------------------------
194+
// - Planet BCs are in here too, can be refactored out
195+
// - Dipole grid must use planet BCs, for now.
196+
// - This kind-of assumes nGCs=2, so may need to be updated.
197+
// - If the first_lower_gc is at iAlt = 1, this may cause issues.
198+
// - The equator-most (j-hat) grid cell will be entirely below min_alt!
199+
for (int iLon = nGCs; iLon < nLons - nGCs; iLon++) {
200+
for (int iLat = nGCs; iLat < nLats - nGCs; iLat++) {
201+
for (iSpecies = 0; iSpecies < nSpecies; iSpecies++) {
187202

188-
if (bcsType == "planet") {
189-
190-
report.print(2, "setting lower bcs to planet");
191-
192-
// Set the lower boundary condition in the last ghost cell:
193-
for (iSpecies = 0; iSpecies < nSpecies; iSpecies++) {
194-
species[iSpecies].density_scgc.slice(nGCs - 1).
195-
fill(species[iSpecies].lower_bc_density);
196-
}
197-
198-
temperature_scgc.slice(nGCs - 1).fill(initial_temperatures[0]);
199-
didWork = true;
200-
}
201-
202-
// fill the second+ grid cells with the bottom temperature:
203-
for (iAlt = nGCs - 2; iAlt >= 0; iAlt--)
204-
temperature_scgc.slice(iAlt) = temperature_scgc.slice(iAlt + 1);
205-
206-
arma_mat sh_ave;
207-
208-
// fill the lower ghost cells with a hydrostatic solution:
209-
for (iSpecies = 0; iSpecies < nSpecies; iSpecies++) {
210-
for (iAlt = nGCs - 2; iAlt >= 0; iAlt--) {
211-
sh_ave =
212-
(species[iSpecies].scale_height_scgc.slice(iAlt) +
213-
species[iSpecies].scale_height_scgc.slice(iAlt + 1)) / 2;
214-
215-
species[iSpecies].density_scgc.slice(iAlt) =
216-
temperature_scgc.slice(iAlt + 1) /
217-
temperature_scgc.slice(iAlt) %
218-
species[iSpecies].density_scgc.slice(iAlt + 1) %
219-
exp(grid.dk_edge_m.slice(iAlt) / sh_ave);
220-
}
221-
222-
for (iAlt = nGCs - 1; iAlt >= 0; iAlt--) {
223-
//std::cout << "before project : " << iAlt << " " << iSpecies << " "
224-
// << species[iSpecies].velocity_vcgc[2](10,10,2) << "\n";
225-
species[iSpecies].velocity_vcgc[2].slice(iAlt) =
226-
species[iSpecies].velocity_vcgc[2].slice(iAlt + 1);
227-
//project_onesided_alt_3rd(species[iSpecies].velocity_vcgc[2], grid, iAlt);
203+
// k-index of 1st lower ghost cell is not constant on the dipole grid.
204+
// On the latlon grid with nGCS=2, this will be 1
205+
iAlt = grid.first_lower_gc(iLon, iLat);
206+
207+
//-----------------------------------------------
208+
// Planet BCs - set to fixed constant values.
209+
//-----------------------------------------------
210+
if (bcsType == "planet" || grid.IsDipole) {
211+
212+
// Fill all lower ghost cells density with lower boundary condition:
213+
species[iSpecies].density_scgc.subcube(iLon, iLat, 0,
214+
iLon, iLat, iAlt - 1).fill(
215+
species[iSpecies].lower_bc_density);
216+
// only fill 1st GC with lower temperature
217+
temperature_scgc(iLon, iLat, iAlt) = initial_temperatures[0];
218+
} // planet bc type
219+
220+
// Set all lower ghost cells to bottom temperature:
221+
temperature_scgc.subcube(iLon, iLat, 0, iLon, iLat, iAlt - 1).fill(
222+
temperature_scgc(iLon, iLat, iAlt));
223+
224+
225+
// 1st ghost cell density is filled with a hydrostatic solution.
226+
sh_ave = (species[iSpecies].scale_height_scgc(iLon, iLat, iAlt)
227+
+ species[iSpecies].scale_height_scgc(iLon, iLat, iAlt + 1)) / 2;
228+
229+
species[iSpecies].density_scgc(iLon, iLat, iAlt) =
230+
temperature_scgc(iLon, iLat, iAlt + 1)
231+
/ temperature_scgc(iLon, iLat, iAlt)
232+
* species[iSpecies].density_scgc(iLon, iLat, iAlt + 1)
233+
* exp(-grid.dk_edge_m(iLon, iLat, iAlt) / sh_ave);
234+
235+
// Vertical velocities: (In GITM this projected down with mesh coeffs)
236+
// Take lowest physical cell's vertical velocity and project it down nGCs cells.
237+
// All "GCs" lower than that have 0 vertical velocity since they're nonphysical.
238+
species[iSpecies].velocity_vcgc[2].subcube(
239+
iLon, iLat, iAlt - 1, size(1, 1, nGCs)).fill(
240+
species[iSpecies].velocity_vcgc[2](iLon, iLat, iAlt + 1));
241+
//project_onesided_alt_3rd(species[iSpecies].velocity_vcgc[2], grid, iAlt);
242+
243+
if (iAlt > nGCs - 1) { // Fill all lower GCs w/ zero vertical velocity
244+
species[iSpecies].velocity_vcgc[2].subcube(iLon, iLat, 0,
245+
size(1, 1, iAlt - 1)).zeros();
246+
247+
}
248+
}
228249
}
229250
}
230251

231-
// Force vertical velocities to be zero in the ghost cells:
232-
for (iDir = 0; iDir < 2; iDir++) {
233-
for (iAlt = 0; iAlt < nGCs; iAlt++) {
234-
for (iSpecies = 0; iSpecies < nSpecies; iSpecies++) {
235-
// species velocity:
236-
species[iSpecies].velocity_vcgc[iDir].slice(iAlt).zeros();
237-
}
252+
didWork = true;
238253

239-
// bulk velocity:
240-
//velocity_vcgc[iDir].slice(iAlt).zeros();
241-
}
242-
}
243254

244255
calc_bulk_velocity();
245256

@@ -281,8 +292,8 @@ bool Neutrals::set_horizontal_bcs(int64_t iDir, Grid grid) {
281292
for (iY = 0; iY < nY; iY++) {
282293
// Constant Gradient for Temperature:
283294
temperature_scgc.tube(iX, iY) =
284-
2 * temperature_scgc.tube(iX - 1, iY) -
285-
temperature_scgc.tube(iX - 2, iY);
295+
2 * temperature_scgc.tube(iX - 1, iY) -
296+
temperature_scgc.tube(iX - 2, iY);
286297

287298
// Constant Value for Velocity:
288299
for (iV = 0; iV < 3; iV++)
@@ -291,8 +302,8 @@ bool Neutrals::set_horizontal_bcs(int64_t iDir, Grid grid) {
291302
// Constant Gradient for densities:
292303
for (int iSpecies = 0; iSpecies < nSpecies; iSpecies++)
293304
species[iSpecies].density_scgc.tube(iX, iY) =
294-
2 * species[iSpecies].density_scgc.tube(iX - 1, iY) -
295-
species[iSpecies].density_scgc.tube(iX - 2, iY);
305+
2 * species[iSpecies].density_scgc.tube(iX - 1, iY) -
306+
species[iSpecies].density_scgc.tube(iX - 2, iY);
296307
}
297308
}
298309
}
@@ -303,8 +314,8 @@ bool Neutrals::set_horizontal_bcs(int64_t iDir, Grid grid) {
303314
for (iY = 0; iY < nY; iY++) {
304315
// Constant Gradient for Temperature:
305316
temperature_scgc.tube(iX, iY) =
306-
2 * temperature_scgc.tube(iX + 1, iY) -
307-
temperature_scgc.tube(iX + 2, iY);
317+
2 * temperature_scgc.tube(iX + 1, iY) -
318+
temperature_scgc.tube(iX + 2, iY);
308319

309320
// Constant Value for Velocity:
310321
for (iV = 0; iV < 3; iV++)
@@ -313,8 +324,8 @@ bool Neutrals::set_horizontal_bcs(int64_t iDir, Grid grid) {
313324
// Constant Gradient for densities:
314325
for (int iSpecies = 0; iSpecies < nSpecies; iSpecies++)
315326
species[iSpecies].density_scgc.tube(iX, iY) =
316-
2 * species[iSpecies].density_scgc.tube(iX + 1, iY) -
317-
species[iSpecies].density_scgc.tube(iX + 2, iY);
327+
2 * species[iSpecies].density_scgc.tube(iX + 1, iY) -
328+
species[iSpecies].density_scgc.tube(iX + 2, iY);
318329
}
319330
}
320331
}
@@ -325,8 +336,8 @@ bool Neutrals::set_horizontal_bcs(int64_t iDir, Grid grid) {
325336
for (iY = nX - nGCs; iY < nY; iY++) {
326337
// Constant Gradient for Temperature:
327338
temperature_scgc.tube(iX, iY) =
328-
2 * temperature_scgc.tube(iX, iY - 1) -
329-
temperature_scgc.tube(iX, iY - 2);
339+
2 * temperature_scgc.tube(iX, iY - 1) -
340+
temperature_scgc.tube(iX, iY - 2);
330341

331342
// Constant Value for Velocity:
332343
for (iV = 0; iV < 3; iV++)
@@ -335,8 +346,8 @@ bool Neutrals::set_horizontal_bcs(int64_t iDir, Grid grid) {
335346
// Constant Gradient for densities:
336347
for (int iSpecies = 0; iSpecies < nSpecies; iSpecies++)
337348
species[iSpecies].density_scgc.tube(iX, iY) =
338-
2 * species[iSpecies].density_scgc.tube(iX, iY - 1) -
339-
species[iSpecies].density_scgc.tube(iX, iY - 2);
349+
2 * species[iSpecies].density_scgc.tube(iX, iY - 1) -
350+
species[iSpecies].density_scgc.tube(iX, iY - 2);
340351
}
341352
}
342353
}
@@ -347,8 +358,8 @@ bool Neutrals::set_horizontal_bcs(int64_t iDir, Grid grid) {
347358
for (iY = nGCs - 1; iY >= 0; iY--) {
348359
// Constant Gradient for Temperature:
349360
temperature_scgc.tube(iX, iY) =
350-
2 * temperature_scgc.tube(iX, iY + 1) -
351-
temperature_scgc.tube(iX, iY + 2);
361+
2 * temperature_scgc.tube(iX, iY + 1) -
362+
temperature_scgc.tube(iX, iY + 2);
352363

353364
// Constant Value for Velocity:
354365
for (iV = 0; iV < 3; iV++)
@@ -357,8 +368,8 @@ bool Neutrals::set_horizontal_bcs(int64_t iDir, Grid grid) {
357368
// Constant Gradient for densities:
358369
for (int iSpecies = 0; iSpecies < nSpecies; iSpecies++)
359370
species[iSpecies].density_scgc.tube(iX, iY) =
360-
2 * species[iSpecies].density_scgc.tube(iX, iY + 1) -
361-
species[iSpecies].density_scgc.tube(iX, iY + 2);
371+
2 * species[iSpecies].density_scgc.tube(iX, iY + 1) -
372+
species[iSpecies].density_scgc.tube(iX, iY + 2);
362373
}
363374
}
364375
}

‎src/neutrals_ics.cpp

Lines changed: 12 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -99,7 +99,7 @@ bool Neutrals::initial_conditions(Grid grid,
9999
std::cout << " NOT Found in MSIS - setting constant\n";
100100

101101
species[iSpecies].density_scgc.slice(0).
102-
fill(species[iSpecies].lower_bc_density);
102+
fill(species[iSpecies].lower_bc_density);
103103
fill_with_hydrostatic(iSpecies, 1, nAlts, grid);
104104
}
105105

@@ -129,8 +129,6 @@ bool Neutrals::initial_conditions(Grid grid,
129129
arma_vec alt1d(nAlts);
130130
arma_vec temp1d(nAlts);
131131

132-
arma_mat H2d(nLons, nLats);
133-
134132
if (nInitial_temps > 0) {
135133
for (iLon = 0; iLon < nLons; iLon++) {
136134
for (iLat = 0; iLat < nLats; iLat++) {
@@ -155,7 +153,7 @@ bool Neutrals::initial_conditions(Grid grid,
155153
iA++;
156154

157155
iA--;
158-
// alt will be between iA and iA+1:
156+
// alt will be between iA and iA+1
159157
r = (alt - initial_altitudes[iA]) /
160158
(initial_altitudes[iA + 1] - initial_altitudes[iA]);
161159
temp1d[iAlt] =
@@ -169,14 +167,18 @@ bool Neutrals::initial_conditions(Grid grid,
169167
}
170168
}
171169
} else
172-
temp1d = 200.0;
170+
temperature_scgc.fill(200.0);
173171

174172
// Make the initial condition in the lower ghost cells to be consistent
175-
// with the actual lowwer BC:
176-
177-
for (int iSpecies = 0; iSpecies < nSpecies; iSpecies++) {
178-
species[iSpecies].density_scgc.slice(0).
179-
fill(species[iSpecies].lower_bc_density);
173+
// with the actual lower BC:
174+
for (iLon = 0; iLon < nLons; iLon ++) {
175+
for (iLat = 0; iLat < nLats; iLat++) {
176+
for (int iSpecies = 0; iSpecies < nSpecies; iSpecies++) {
177+
species[iSpecies].density_scgc.subcube(
178+
iLon, iLat, 0, iLon, iLat, grid.first_lower_gc(iLon, iLat)).fill(
179+
species[iSpecies].lower_bc_density);
180+
}
181+
}
180182
}
181183

182184
report.print(2, "Calculating scale height");

0 commit comments

Comments
 (0)