Skip to content

Conversation

@cscjlan
Copy link

@cscjlan cscjlan commented Feb 17, 2025

FsGrid

Vlasiator has been updated to use the refactored FsGrid. This is visible everywhere FsGrid was used in. Some notable differences are dicussed below.

First, there's only a single FsGrid object and it's called fsgrid. It no longer stores any data. The data that used to be stored by the different grids are now owned by a structure called FsData, one for each ex-grid. The data is not passed around by passing the FsData structures, but rather by passing an std::span. You can get an std::span from FsData with the view() method: auto technicalSpan = technicalData.view();

updateGhostCells now takes in an std::span: fsgrid.updateGhostCells(technical); instead of technicalGrid.updateGhostCells().

Data is no longer accessed through the grid.get(i, j, k) method, but rather by contructing an FsStencil, and using it to compute a 1D index from a 3D index:

auto stencil = fsgrid.makeStencil(i, j, k);
// Get the center element ('ooo'), i.e. the one corresponding to i, j, k
const auto center = technical[stencil.ooo()];

// The element in the -x direction from the center, i.e. i - 1, j, k
const auto x_left = technical[stencil.moo()];

This has the benefit that some values that are constant throughout the simulation can be computed once and stored in the fsgrid object. A single FsStencil is constructed for each iteration of a loop (i.e. one for each i, j, k). The constant values can be accessed by the FsStencil when it computes the 1D value from the 3D index it was constructed with.

Fieldsolver

Fieldsolver has been refactored to use the updated FsGrid. Some files have undergone heavier refactoring (ldz_electric_field.cpp), some less heavy. The underlying reason for all of the refactoring is to make it easier to change the fieldsolver to use the updated FsGrid: with less code repetition there's less work to be done (and thus fewer bugs) when changing the code to use the stencil and the spans. An additional benefit was a speed up especially in the electric field computation.

derivatives.cpp and ldz_volume.cpp use the newly added fsgrid.parallel_for, which is an architecture/processor independent way of doing the desired computation. It takes in a lambda, which it calls for each local cell on the fsgrid.

… of dmoments::dPed was computed regardless of shouldCalculateMoments flag
@ykempf
Copy link
Contributor

ykempf commented Aug 12, 2025

All right I seem to have botched something(s). However nominally the last 3 commits should have dealt with

  • getting rid of the numerous variants of parallel_for_* loop interfaces
  • removing RNG inits for perturbed B
  • removing capturing gridSpacing by using the newly passed coordinates.

@ykempf
Copy link
Contributor

ykempf commented Aug 12, 2025

Oh yeah forgot to push out the fsgrid commit. 🫠

@ykempf
Copy link
Contributor

ykempf commented Aug 20, 2025

The rabbit hole of #1161 convinced me that the testpackage diffs are very unlikely to have arisen from this PR. I will thus proceed with the remaining open things here.

Copy link
Contributor

@ursg ursg left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Incomplete review so far, I'll go through some more later. But functional approval, of course.
I'm not fully happy with the style in a number of places, but those are things that can... settle over time.

Comment on lines +132 to +136
fsgrid.parallel_for([](int timerId) -> phiprof::Timer { return phiprof::Timer{timerId}; },
phiprof::initializeTimer("setBackgroundFieldToZero"), technical,
[=](const fsgrid::Coordinates &coordinates, const fsgrid::FsStencil& stencil, cuint sysBoundaryFlag, cuint sysBoundaryLayer) {
for (size_t i = 0; i < bgb[stencil.ooo()].size(); i++) {
bgb[stencil.ooo()][i] = 0.0;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This parallel for has inconsistent indentation with the others above.

Comment on lines +51 to +52
if (P::systemWriteAllDROs || lowercase == "fg_b" ||
lowercase == "b") { // Bulk magnetic field at Yee-Lattice locations
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm completely ok with most of your autoreformattings in this PR, but breaking all these lines in the datareducer file actually makes the whole source file harder to read. Is this something you can semi-easily revert by hand, or should I do that?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These were all done automatically by Juhana's software, according to the repo's style file. In things I did manually, I tried to stick to some sort of "what was done so far in the team", but I did not revert Juhana's automatic formatting. I can provide the style file I came up with when trying to improve it. I would really strongly suggest from this side of the bay that pliiis, deploy some CI tool doing all formatting for everyone, or at least a style file that works for all. However I stumbled upon options that changed in recentish versions of clang formatting, and/or git-clang-format was not available on major platforms we use (LUMI), so I think the more robust way would be a CI pre-commit or similar feature holding the ground truth and making everything consistent. At the price of having to create that ground truth file...

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

As @ykempf mentioned in a comment down the line, many of these could be "reverted" quite easily by increasing the allowed line width.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

"They're more like guidelines anyway".

I'm all in favour of allowing for longer lines here, because the readability actually benefits from the better formatting in many cases.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So! Set up some way or other for this to be rolled out automatically, fix the formatting options to your desire, and off we go! Some options allow or disallow e.g. one parameter per line when breaking, but also some options were recent enough not to be available on all platforms, so I suppose having one central reference instance will be easiest to maintain (I tried installing git-clang-format on some machine without having to build the whole of LLVM and failed with that...).

Comment on lines +250 to +253
outputReducer->addOperator(new DRO::DataReductionOperatorPopulations<Real>(
pop + "/vg_rho", i, offsetof(spatial_cell::Population, RHO), 1));
outputReducer->addMetadata(outputReducer->size() - 1, "1/m^3", "$\\mathrm{m}^{-3}$",
"$n_\\mathrm{" + pop + "}$", "1.0");
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

likewise, breaking these unit lines just makes things unneccesarily hard to read.

Comment on lines +874 to +909
outputReducer->addOperator(new DRO::DataReductionOperatorBVOLDerivatives("vg_derivatives/vg_dperbxvoldx",
bvolderivatives::dPERBXVOLdx, 1));
outputReducer->addOperator(new DRO::DataReductionOperatorBVOLDerivatives("vg_derivatives/vg_dperbxvoldy",
bvolderivatives::dPERBXVOLdy, 1));
outputReducer->addOperator(new DRO::DataReductionOperatorBVOLDerivatives("vg_derivatives/vg_dperbxvoldz",
bvolderivatives::dPERBXVOLdz, 1));
outputReducer->addOperator(new DRO::DataReductionOperatorBVOLDerivatives("vg_derivatives/vg_dperbyvoldx",
bvolderivatives::dPERBYVOLdx, 1));
outputReducer->addOperator(new DRO::DataReductionOperatorBVOLDerivatives("vg_derivatives/vg_dperbyvoldy",
bvolderivatives::dPERBYVOLdy, 1));
outputReducer->addOperator(new DRO::DataReductionOperatorBVOLDerivatives("vg_derivatives/vg_dperbyvoldz",
bvolderivatives::dPERBYVOLdz, 1));
outputReducer->addOperator(new DRO::DataReductionOperatorBVOLDerivatives("vg_derivatives/vg_dperbzvoldx",
bvolderivatives::dPERBZVOLdx, 1));
outputReducer->addOperator(new DRO::DataReductionOperatorBVOLDerivatives("vg_derivatives/vg_dperbzvoldy",
bvolderivatives::dPERBZVOLdy, 1));
outputReducer->addOperator(new DRO::DataReductionOperatorBVOLDerivatives("vg_derivatives/vg_dperbzvoldz",
bvolderivatives::dPERBZVOLdz, 1));
outputReducer->addMetadata(outputReducer->size() - 9, "T/m", "$\\mathrm{T}\\,\\mathrm{m}^{-1}$",
"$\\Delta B_{X,\\mathrm{per,vol,vg}} (\\Delta X)^{-1}$", "1.0");
outputReducer->addMetadata(outputReducer->size() - 8, "T/m", "$\\mathrm{T}\\,\\mathrm{m}^{-1}$",
"$\\Delta B_{X,\\mathrm{per,vol,vg}} (\\Delta Y)^{-1}$", "1.0");
outputReducer->addMetadata(outputReducer->size() - 7, "T/m", "$\\mathrm{T}\\,\\mathrm{m}^{-1}$",
"$\\Delta B_{X,\\mathrm{per,vol,vg}} (\\Delta Z)^{-1}$", "1.0");
outputReducer->addMetadata(outputReducer->size() - 6, "T/m", "$\\mathrm{T}\\,\\mathrm{m}^{-1}$",
"$\\Delta B_{Y,\\mathrm{per,vol,vg}} (\\Delta X)^{-1}$", "1.0");
outputReducer->addMetadata(outputReducer->size() - 5, "T/m", "$\\mathrm{T}\\,\\mathrm{m}^{-1}$",
"$\\Delta B_{Y,\\mathrm{per,vol,vg}} (\\Delta Y)^{-1}$", "1.0");
outputReducer->addMetadata(outputReducer->size() - 4, "T/m", "$\\mathrm{T}\\,\\mathrm{m}^{-1}$",
"$\\Delta B_{Y,\\mathrm{per,vol,vg}} (\\Delta Z)^{-1}$", "1.0");
outputReducer->addMetadata(outputReducer->size() - 3, "T/m", "$\\mathrm{T}\\,\\mathrm{m}^{-1}$",
"$\\Delta B_{Z,\\mathrm{per,vol,vg}} (\\Delta X)^{-1}$", "1.0");
outputReducer->addMetadata(outputReducer->size() - 2, "T/m", "$\\mathrm{T}\\,\\mathrm{m}^{-1}$",
"$\\Delta B_{Z,\\mathrm{per,vol,vg}} (\\Delta Y)^{-1}$", "1.0");
outputReducer->addMetadata(outputReducer->size() - 1, "T/m", "$\\mathrm{T}\\,\\mathrm{m}^{-1}$",
"$\\Delta B_{Z,\\mathrm{per,vol,vg}} (\\Delta Z)^{-1}$", "1.0");
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also this is unnecessary, IMHO

Comment on lines +2286 to +2288
Real theta = acos(grid.nodes[n].x[2] / sqrt(grid.nodes[n].x[0] * grid.nodes[n].x[0] +
grid.nodes[n].x[1] * grid.nodes[n].x[1] +
grid.nodes[n].x[2] * grid.nodes[n].x[2])); // Latitude
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Here's an exception, that I actually appreciate! :)

Comment on lines +2614 to +2619
Real deltaE =
(particle_energy[p + 1] - particle_energy[p]) * 1e3 * physicalconstants::CHARGE; // dE in J
numberFlux += grid.nodes[i].parameters[ionosphereParameters::RHON] *
sqrt(1. / (2. * M_PI * physicalconstants::MASS_ELECTRON)) * particle_energy[p] /
temp_keV / sqrt(temp_keV * 1e3 * physicalconstants::CHARGE) * deltaE *
exp(-energyparam); // Flux 1/m^2/s
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same here.

Comment on lines -812 to +586
dPsq = std::max((pow(myP[0] - otherP[0], 2) + pow(myP[1] - otherP[1], 2) + pow(myP[2] - otherP[2], 2)) / (2 * myRho * maxU), dPsq);
dPsq = std::max((pow(myP[0] - otherP[0], 2) + pow(myP[1] - otherP[1], 2) + pow(myP[2] - otherP[2], 2)) /
(2 * myRho * maxU),
dPsq);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Here's another particularly unpretty formatting example.

Comment on lines +112 to +140
perturbedResult[Rec::a_yy] =
HALF * (der_i2j1k1[fsgrids::dperb::dPERBxdyy] + der_i1j1k1[fsgrids::dperb::dPERBxdyy]);
perturbedResult[Rec::a_zz] =
HALF * (der_i2j1k1[fsgrids::dperb::dPERBxdzz] + der_i1j1k1[fsgrids::dperb::dPERBxdzz]);
perturbedResult[Rec::a_yz] =
HALF * (der_i2j1k1[fsgrids::dperb::dPERBxdyz] + der_i1j1k1[fsgrids::dperb::dPERBxdyz]);
perturbedResult[Rec::a_xyy] = (der_i2j1k1[fsgrids::dperb::dPERBxdyy] - der_i1j1k1[fsgrids::dperb::dPERBxdyy]);
perturbedResult[Rec::a_xyz] = (der_i2j1k1[fsgrids::dperb::dPERBxdyz] - der_i1j1k1[fsgrids::dperb::dPERBxdyz]);
perturbedResult[Rec::a_xzz] = (der_i2j1k1[fsgrids::dperb::dPERBxdzz] - der_i1j1k1[fsgrids::dperb::dPERBxdzz]);

perturbedResult[Rec::b_xx] =
HALF * (der_i1j2k1[fsgrids::dperb::dPERBydxx] + der_i1j1k1[fsgrids::dperb::dPERBydxx]);
perturbedResult[Rec::b_xz] =
HALF * (der_i1j2k1[fsgrids::dperb::dPERBydxz] + der_i1j1k1[fsgrids::dperb::dPERBydxz]);
perturbedResult[Rec::b_zz] =
HALF * (der_i1j2k1[fsgrids::dperb::dPERBydzz] + der_i1j1k1[fsgrids::dperb::dPERBydzz]);
perturbedResult[Rec::b_xxy] = (der_i1j2k1[fsgrids::dperb::dPERBydxx] - der_i1j1k1[fsgrids::dperb::dPERBydxx]);
perturbedResult[Rec::b_xyz] = (der_i1j2k1[fsgrids::dperb::dPERBydxz] - der_i1j1k1[fsgrids::dperb::dPERBydxz]);
perturbedResult[Rec::b_yzz] = (der_i1j2k1[fsgrids::dperb::dPERBydzz] - der_i1j1k1[fsgrids::dperb::dPERBydzz]);

perturbedResult[Rec::c_xx] =
HALF * (der_i1j1k2[fsgrids::dperb::dPERBzdxx] + der_i1j1k1[fsgrids::dperb::dPERBzdxx]);
perturbedResult[Rec::c_xy] =
HALF * (der_i1j1k2[fsgrids::dperb::dPERBzdxy] + der_i1j1k1[fsgrids::dperb::dPERBzdxy]);
perturbedResult[Rec::c_yy] =
HALF * (der_i1j1k2[fsgrids::dperb::dPERBzdyy] + der_i1j1k1[fsgrids::dperb::dPERBzdyy]);
perturbedResult[Rec::c_xxz] = (der_i1j1k2[fsgrids::dperb::dPERBzdxx] - der_i1j1k1[fsgrids::dperb::dPERBzdxx]);
perturbedResult[Rec::c_xyz] = (der_i1j1k2[fsgrids::dperb::dPERBzdxy] - der_i1j1k1[fsgrids::dperb::dPERBzdxy]);
perturbedResult[Rec::c_yyz] = (der_i1j1k2[fsgrids::dperb::dPERBzdyy] - der_i1j1k1[fsgrids::dperb::dPERBzdyy]);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

And here, only ~50% of the lines have been broken now.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

(This is anyway code that probably nobody will ever read through with any intent or detail... so why isn't this just autogenerated anyway?)

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is inherited directly from Arto's hand-crafted first implementation out of the paper, I/we only changed the arrays, namespaces etc in hiostory but the indices are still as typed by Him in the olden days.

Comment on lines +65 to +79
bool propagateFields(fsgrids::perbspan perb,
fsgrids::perbspan perbdt2,
fsgrids::efieldspan e,
fsgrids::efieldspan edt2,
fsgrids::ehallspan ehall,
fsgrids::egradpespan egradpe,
fsgrids::egradpespan egradpedt2,
fsgrids::momentsspan moments,
fsgrids::momentsspan momentsdt2,
fsgrids::dperbspan dperb,
fsgrids::dmomentsspan dmoments,
fsgrids::dmomentsspan dmomentsdt2,
fsgrids::bgbspan bgb,
fsgrids::volspan vol,
fsgrids::technicalspan technical, FieldSolverGrid &fsgrid, SysBoundary& sysBoundaries,
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why aren't all of these spans just collected together into one proper "fsgrids" object by now? I mean, it's great that they live in their own namespace, but for all intents and purposes, they might as well be handed back and forth in a single pointer, if they were sharing a common struct.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

(Also, the fieldsolver functions could then be methods of that object)

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

And now I saw further down that FsGridSolverData does exist. So how about passing that one here?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think so far the idea was to pass the needed spans down the chain, which can transparently point to device or host memory in the end.

Copy link
Author

@cscjlan cscjlan Sep 15, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

IIRC I kept the function signatures on the vlasiator side as close to the original as possible if it wasn't necessary to change them. But definitely one could just pass FieldSolverData here.


I would keep the functionality separate from the FieldSolverData object, i.e. not implement any methods for it. First of all, the FieldSolverData doesn't even own the data, it's just a wrapper over multiple std::spans -- i.e. pointers & lengths -- the original purpose of which was to make the tens of lambda expressions simpler in datareducer.cpp to simplifiy changing the underlying objects, because they can be changed in one place.

Secondly keeping data and functionality separated was one of the major points of the refactoring on the FsGrid side. This way e.g. the initialization of FsData and FieldSolverGrid can happen separately from each other and there are no partially initialized objects: "on this line the MPI initialization has happened for object A, but the data is garbage, but after these lines also the data has been initialized for object A". With separation, FieldSolverGrid is a fully functioning object after construction and likewise for each FsData (albeit they're full of zeroes).

This separation also simplifies e.g. testing of FsData: when it's implemented to work on the GPU, testing it can be done without any references to MPI/OpenMP, since it neither knows nor cares about those.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I see. Ok then, given what we have learned about memory handling in hybrid architectures, it makes a lot of sense to handle the spans individually and carefully.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hooray, my least favourite source file got even messier now. :)

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Show it some love, come on! I assume the hyperresistivity one will not be much different. :)

@ykempf
Copy link
Contributor

ykempf commented Sep 11, 2025

By my eye, most if not all of the lines you disliked in terms of formatting are due to the max character width setting in the config file.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

5 participants