ammarhakim / gkylzero

Lowest, compiled layer of Gkeyll infrastructure.
MIT License
22 stars 5 forks source link

[DR] Bringing Vlasov Poisson into the Vlasov app #456

Closed manauref closed 3 weeks ago

manauref commented 1 month ago

Presently we have a separate Vlasov Poisson app in the gk-g0-app_vp branch (and some sub branches). We'd like to bring this into the main Vlasov app for (at least) 2 reasons:

  1. So it's easier to use components common to both solvers (e.g. collisions).
  2. So that VP undergoes the same level of scrutiny as other solvers, giving users more peace of mind.

The approach to do so is as follows:

  1. A user will indicate in the app input table (gkyl_vm) that they wish to run an electrostatic simulation with is_electrostatic, such as

    struct gkyl_vm vm = {
    .name = "vp_sheath_1x1v_p2",
    
    .cdim = ctx.cdim, .vdim = ctx.vdim,
    .lower = { ctx.x_min },
    .upper = { ctx.x_max },
    .cells = { ctx.Nx },
    .poly_order = ctx.poly_order,
    .basis_type = app_args.basis_type,
    
    .is_electrostatic = true,
    ...
  2. Additionally, input files will specify alternative entries in the field input table (gkyl_vlasov_field) like
    struct gkyl_vlasov_field field = {
    .epsilon0 = ctx.epsilon0,
    .poisson_bcs = {
      .lo_type = { GKYL_POISSON_DIRICHLET },
      .up_type = { GKYL_POISSON_NEUMANN },
      .lo_value = { 0.0 }, .up_value = { 0.0 }
    },
    // External potentials are optional:
    .ext_potentials = ext_potentials_profiles,
    .ext_potentials_ctx = &ctx,,
    .ext_potentials_evolve = false,
    };
  3. Internally, the only new object we will create is a apps/vp_field.c, which handles the Poisson problem and the external potentials. The Vlasov app will choose between vm_field or vp_field based on the is_electrostatic input (default=false).
  4. In order to not further balloon the Vlasov Maxwell solver, to avoid its accidental modification, and make the code more modular, we will create a separate rk3 stepping function for VP. This update function will be chosen based on is_electrostatic in a similar way to how we choose the operator split update (e.g. for implicit BGK), e.g.
    if (app->is_electrostatic) {
    if (app->has_implicit_coll_scheme || app->has_fluid_em_coupling) {
      app->update_func = vlasov_poisson_update_op_split;
    }
    else {
      app->update_func = vlasov_poisson_update_ssp_rk3;
    }
    }
    else {
    if (app->has_implicit_coll_scheme || app->has_fluid_em_coupling) {
      app->update_func = vlasov_update_op_split;
    }
    else {
      app->update_func = vlasov_update_ssp_rk3;
    }
    }

    The reason for this is that VP uses a different stepping logic than VM, and is more like that used in GK. A version of vlasov_poisson_update_ssp_rk3 already exists and has been validated in the gk-g0-app_vp branch.

  5. The final app-related point, which has not been fully worked out, pertains field diagnostics. Certain field diagnostics like the field energy diagnostic needs to be modified. I think on the user-side we'll keep the call
    gkyl_vlasov_app_calc_field_energy(gkyl_vlasov_app* app, double tm)

    which internally calls

    vm_field_calc_energy(app, tm, app->field);

    So what we'll do is change the user-facing function so it actually does

    app->field_calc_energy(app, tm, app->field);

    and this pointer will be set to vm_field_calc_energy or vp_field_calc_energy based on is_electrostatic.

  6. Lastly, we will separate VP from the underlying vlasov DG solver (dg_updater_vlasov), creating dg_updater_vlasov_poisson, so that its kernels and collisionless solvers are isolated for future modification without tampering with the VM solver (and its various flavors).

Much of this work is prototyped in the gk-g0-app_vp branch.

manauref commented 3 weeks ago

Done in PR #473