Sub-grid model for pitch-angle diffusion #846
Conversation
Remove abs value for diff coeff.
Added include<string> in various files Built Spanier diffusion Moved old diffusion
…fixed empty cells.
|
Hmm, CI compilation fails: We've so far always used version 1 of vectorclass, will moving to version 2 cause some new issues? |
| for (uint k = 0; k < WID; ++k) for (uint j = 0; j < WID; ++j) { // Iterate through coordinates (z,y) | ||
|
|
||
| //Get velocity space coordinates | ||
| const Vec4d VX(parameters[n * BlockParams::N_VELOCITY_BLOCK_PARAMS + BlockParams::VXCRD] |
There was a problem hiding this comment.
This is now hard-coded to vec4d. It's a velocity, not a distribution function value, so it should be Real instead of Realf, but hard-coding the length to 4 will be messy in the future. I would recommend instead doing something like looping j over vec_per_plane - look in acceleration for suggestions on how to accomplish this.
There was a problem hiding this comment.
You can also use #pragma omp simd here to suggest vectorisation on the compiler level (and run some tests if it works or not)
There was a problem hiding this comment.
We attempted #pragma omp simd autovectorization here, but the compiler is adamant that a memory access prevents vectorization altogether:
vlasovsolver/velocity_space_diffusion.cpp:241:48: missed: couldn't vectorize loop
vlasovsolver/velocity_space_diffusion.cpp:276:29: missed: not vectorized: no vectype for stmt: CellValue_803 = *_802;
scalar_type: Realf
I suppose there would be a way to nudge the compiler in the right direction here, but this probably needs a special kind of witchcraft that I am not familiar with.
| {cell.parameters[CellParams::P_12],cell.parameters[CellParams::P_22],cell.parameters[CellParams::P_23]}, | ||
| {cell.parameters[CellParams::P_13],cell.parameters[CellParams::P_23],cell.parameters[CellParams::P_33]}}; | ||
|
|
||
| // Build Rotation Matrix |
There was a problem hiding this comment.
Yeah this looks pretty horrible :D
it would be a lot better to use the Eigen library for this.
Replace muspace access with a horribly convoluted C language construct.
| @@ -227,7 +227,12 @@ void computeNewTimeStep(dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>& mpi | |||
| Real subcycleDt; | |||
|
|
|||
| // reduce/increase dt if it is too high for any of the three propagators or too low for all propagators | |||
| if ((P::dt > dtMaxGlobal[0] * P::vlasovSolverMaxCFL || | |||
| if (P::dt == P::dt_ceil) { | |||
| isChanged = false; | |||
There was a problem hiding this comment.
Note: This means that if the timestep ever reaches the ceiling, it can never become smaller again. That's dangerous!
Instead, add ((P::dt_ceil > 0)&&(P::dt>P::dt_ceil)) to the if-clause-list on line 235-240.
There was a problem hiding this comment.
ah now I see the continue-mechanism Yann talked about. Not so straightforward, then...
There was a problem hiding this comment.
Went through the logic again, and in fact the above suggested change should do the trick, if I'm correct.
Enjoy.