Ionization in DSMC

[ancientflow]

I am attempting to simulate the ignition ionization of Xe using DSMC, and I have set the ndt of DSMC to 10. This setting can meet the requirement that the collision probability is much less than 1. However, the electron density result I obtained is much smaller than that given in the literature. By checking the source code, I found that when the SplitAndScatterFunc class is called, the dt parameter does not seem to be input. Does this function currently not support performing collisions at intervals of a certain time step, that is, can ndt only be 1 when processing ionization in DSMC?

[roelof-groenewald]

Hi @ancientflow. The ndt parameter is used at this line:

warpx/Source/Particles/Collision/CollisionBase.cpp

Line 23 in 09637e1

pp_collision_name, "ndt", ndt);

Do you see a change in results when you use ndt = 1?

[ancientflow]

@roelof-groenewald Hi, I have used ndt=1 and ndt=10 to get results. They are almost the same and my problem might lie elsewhere.

[ancientflow]

@roelof-groenewald In this simulation, a spatially distributed boundary condition is needed. I've noticed that this might not be supported at present, so I revised the electrostatic solver code at warpx/Source/FieldSolver/ElectrostaticSolvers/ElectrostaticSolver.cpp line 120. I added the following code:

    static const double ncell = 128, halfncell = ncell / 2;
    static const double gap = 0.005 / ncell,
                        index_sq_min = (0.0021 / 2 / gap) * (0.0021 / 2 / gap),
                        index_sq_max = (0.0031 / 2 / gap) * (0.0031 / 2 / gap);
    auto phi_field = warpx.m_fields.get(FieldType::phi_fp, 0);
    for (MFIter mfi(*phi_field, TilingIfNotGPU()); mfi.isValid(); ++mfi) {
        const Box& box = mfi.tilebox();

        Array4<Real> const& phi = phi_field->array(mfi);
        ParallelFor(box, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
            double i_ex = i - halfncell, j_ex = j - halfncell;
            if (k == domain.smallEnd(2)) {
                if (i_ex * i_ex + j_ex * j_ex >= index_sq_min &&
                    i_ex * i_ex + j_ex * j_ex <= index_sq_max)
                    phi(i, j, k) = 300.0;
                else
                    phi(i, j, k) = 0.0;
            }
        });
    }

The electric potential in the other areas and boundaries is zero. However, when I opened the diagnostic file, I got the result that the area where the potential is set at 300V only had a potential of only approximately 285V. How can I solve this problem?

[roelof-groenewald]

When the fields are output during diagnostics all fields are interpolated to cell centers. During the simulation phi is defined on the cell nodes, so I expect you are seeing slightly lower values since it is the average of the boundary value and the potential on the next grid point inside the domain. Does that make sense?

[ancientflow]

Thanks for your reply. I have confirmed that these are the interpolated values and the boundary values are what I need.