Vlasov equation integration issue (hyper.deal)

[Alexander Kiselyov]

Dear all,

I'm having a mysterious issue with hyper.deal, but I suspect the reason
might lie in my mistakes in using or misunderstanding deal.II.

I'm currently trying to develop a 6-D Vlasov FEM solver using hyper.deal
accounting for both electric and magnetic field as a first step towards
a self-consistent Vlasov-Maxwell solver. The current physical problem is
collisionless plasma in a constant magnetic field B. From my
understanding, spherically symmetric velocity distribution should not
change in time in this case (magnetic field would induce its rotation
around B vector).

Diagnostics tools are density ($\int f(r,v) dv$), velocity distribution
($\int f(r,v) dr$), full f(r,v) output (r, v and f(r,v) at each
quadrature point) and advection vector (transport velocity) field. The
issue is that when using a function depending on v as an advection
vector in the df/dv term of the Vlasov equation artifacts appear (below
the examples are for the Lorenz force without electric field, but it is
observable for something artificial like $(v \cdot B) B$ too). They
manifest themselves as ever increasing values "along" magnetic field for
odd quadrature points in velocity distribution plots (see attached
images). What the strangest thing is that saving solution values at
quadrature points, doing naïve integration (i.e. just summing f values)
and plotting using Gnuplot does _not_ reveal such artifacts: the
velocity distribution remains spherically symmetric. It is puzzling that
their amplitude does _not_ depend on the usual numerical parameters: dt
and grid density.

See `artifact_vtu.png` for velocity distribution calculated using
deal.II/hyper.deal routines, `artifact_naive.png` for a naïve velocity
distribution calculation. Advection vector (transport velocity) fields
for coordinate and velocity spaces are shown in `advection_x.png` at v =
(0.014088, -0.38909, -0.36091) and in `advection_v.png` at r =
(-0.014088, -0.014088, -0.0625). Numerical parameters are dt = 0.01, T =
0.1, the grid is a hypercube [-0.5; 0.5]^6 with the size of 8
subdivisions along each dimension. Magnetic field is constant B = (1, 0,
0), normalized in such a way that $F_{Lorenz} = q [v \times B]$.

The velocity distribution routine is implemented along the same lines as
`hyperdeal::VectorTools::velocity_space_integration`
(`include/vector_tools.h`):

data.template cell_loop<Vector_Out, Vector_In>(
[&](const auto&, Vector_Out& dest, const Vector_In& src_,
const auto cell) mutable
{
phi.reinit(cell);
phi.read_dof_values(src_);
phi_x.reinit(cell.x);
phi_v.reinit(cell.v);

const VectorizedArrayType* data_ptr_src = phi.get_data_ptr();
VectorizedArrayTypeV* data_ptr_dst = phi_v.begin_dof_values();

const auto nv_factor = dealii::Utilities::pow(n_points, dim_v);
const auto N_lanes = phi.n_vectorization_lanes_filled();

for (unsigned int qv = 0; qv < phi_v.n_q_points; ++qv) {
VectorizedArrayTypeV sum_x = 0.0;

for (unsigned int qx = 0; qx < phi_x.n_q_points; ++qx) {
const auto q = qx + qv*nv_factor;
for (unsigned int lane = 0; lane < N_lanes; ++lane)
sum_x += data_ptr_src[q][lane] * phi_x.JxW(qx)[lane];
}

data_ptr_dst[qv] = sum_x;
}

phi_v.distribute_local_to_global(dest);
},
dst,
src
);

The extraction of solution values at quadrature points is quite similar
(`include/output.h`, line 427):

matrix_free.template cell_loop<int, VectorType>(
[&out, &phi, &phi_x, &phi_v](const auto&, int&,
const VectorType& src,
const auto cell) mutable
{
phi.reinit(cell);
phi.read_dof_values(src);
phi_x.reinit(cell.x);
const unsigned int index_v = cell.v / VectorizedArrayTypeV::size();
const unsigned int lane_v = cell.v % VectorizedArrayTypeV::size();
phi_v.reinit(cell.v);

const VectorizedArrayType* data_ptr_src = phi.get_data_ptr();

const auto nv_factor = dealii::Utilities::pow(n_points, dim_v);
const auto N_lanes = phi.n_vectorization_lanes_filled();
for (unsigned int qv = 0; qv < phi_v.n_q_points; ++qv) {
for (unsigned int qx = 0; qx < phi_x.n_q_points; ++qx) {
const auto qx_point = phi_x.quadrature_point(qx);
const auto qv_point = phi_v.quadrature_point(qv);
const auto q_value = data_ptr_src[qx + qv*nv_factor];

for (unsigned int lane = 0; lane < N_lanes; ++lane) {
double cache = 0.0;
// Write coordinates
for (unsigned int d = 0; d < dim_x; ++d) {
// all lanes filled
cache = qx_point[d][lane];
out.write(reinterpret_cast<const char*>(&cache), sizeof cache);
}
for (unsigned int d = 0; d < dim_v; ++d) {
// only one lane used for v
cache = qv_point[d][lane_v];
out.write(reinterpret_cast<const char*>(&cache), sizeof cache);
}
// Write value
cache = q_value[lane];
out.write(reinterpret_cast<const char*>(&cache), sizeof cache);
}
}
}
},
dummy, solution
);

Attached below is also a tarball `plasma-model.tar.bz2` with the
complete buildable C++17 sources. The plots can be recreated by running
`vlasov` binary with `--output-full` argument via `mpirun` (*the output
is about 36 GiB large!*). VTUs can be viewed using `vlasov_state.pvsm`
in ParaView, naïve integration by running `reader -i 10 --integrate` and
then plotting `density.gpi` with Gnuplot.

What advice can you give on debugging the issue? Are there flaws in my
understanding of deal.II/hyper.deal or numerical peculiarities of such a
model?

Best regards,
Alexander Kiselyov

P.S. System specifications:

* A single machine with Intel Core i7-12700K CPU
* NixOS 21.11 amd64 (Linux kernel 5.16.3, gcc 10.3.0, OpenMPI 4.1.1)
* Deal.II 10.0.0-prerelease (b01ae8a)
* Hyper.deal (9bfd4b2)
* ParaView 5.8.0

P.P.S. Also I would like to thank the developers of deal.II and
hyper.deal for making such great projects available as free software.

[Peter Munch]

Hi Alexander,

Nice that you found hyper.deal!

Since this seems to be more a hyper.deal related issue, I would suggest to move the discussion to: https://github.com/hyperdeal/hyperdeal/discussions. 

On the first glance, the code seems to be correct but I have one or two ideas how to simplify the problem so that it is easier to debug!

Thanks,

Peter

[Alexander Kiselyov]

Thank you very much for the reply! Copied my letter to
https://github.com/hyperdeal/hyperdeal/discussions/90.

Regards,
Alexander

> --
> The deal.II project is located at http://www.dealii.org/
> <http://www.dealii.org/>
> For mailing list/forum options, see
> https://groups.google.com/d/forum/dealii?hl=en
> <https://groups.google.com/d/forum/dealii?hl=en>
> ---
> You received this message because you are subscribed to the Google
> Groups "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send
> an email to dealii+un...@googlegroups.com
> <mailto:dealii+un...@googlegroups.com>.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/dealii/98ea374b-848b-491b-8675-8e09fbb08476n%40googlegroups.com
> <https://groups.google.com/d/msgid/dealii/98ea374b-848b-491b-8675-8e09fbb08476n%40googlegroups.com?utm_medium=email&utm_source=footer>.