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 On Thursday, 3 February 2022 at 12:12:20 UTC+1 Alexander Kiselyov wrote: > 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. -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see 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 [email protected]. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/98ea374b-848b-491b-8675-8e09fbb08476n%40googlegroups.com.
