ammarhakim / gkylzero

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

[DR] Ability to specify biased planes in Poisson problems #535

Open manauref opened 4 hours ago

manauref commented 4 hours ago

We wish to rigorously enforce the twist-shift-icity of the potential phi in GK clopen simulations, see issue #504 . The only way we now to do this now (Phases 1-4) all modify the potential after the Poisson solver at the boundary, which introduces a discontinuity in the x direction at z skin cells abutting the LCFS.

@gwhammett pointed out that the only way we can think of avoiding this without messing with the continuity of phi is to enforce a condition phi(x=LCFX, z=zmin,zmax)=constant. This way the TS BC has no effect on it, so enforcing TS(phi) doesn't change phi(x=LCFX, z=zmin,zmax), and phi remains continuous in x.

Here's a proposal for how to do this based on work @Antoinehoff and I did this week:

  1. Create a new struct to specify a biased plane in the Poisson solver, such as
    struct gkyl_poisson_bias_plane {
    int dir; // Direction perpendicular to the plane.
    double loc; // Location of the plane in the 'dir' dimension.
    double val; // Biasing value.
    };

    The user can specify a list of biased planes with

    struct gkyl_poisson_bias_plane_list {
    int num_bias_plane;
    struct gkyl_poisson_bias_plane *bp;
    };

    Note that here "planes" means a point in 1D, a line in 2D, or a surface in 3D.

  2. Implement a new kernel which replaces the <i,j,poisson_value> entry in the triples object, used to populate the LHS matrix of the Poisson equation, corresponding to a node on this plane with <i,j,1> or <i,j,0>, so that the equation for that node becomes phi_k = bias_value (as we do for Dirichlet BCs at the boundary).
  3. Implement another kernel that does the same thing for the RHS source vector of the Poisson problem, but replaces the corresponding row of the vector with bias_value (as we do for Dirichlet BCs at the boundary).
  4. In GK we need to solve 2D problems at planes, and for clopen sims we only need to bias at the lower/upper plane of the domain. We will pass a gkyl_poisson_bias_plane_list to deflated_fem_poisson, but this object will only give it to the lower/upper Poisson solves.

Note about steps 2 & 3: We could consider enforcing these conditions in the same kernels that discretize the Poisson equation, as we do for Dirichlet BCs. However for now I think it's safer to keep this biasing operation in separate code as to not jeopardize the integrity of the Poisson solver we know works.

The fem_poisson changes are being prototyped in gk-g0-app_FLR because there are other fem_poisson changes there we need to build this on top of.

Antoinehoff commented 3 hours ago

Here is an illustration of what the matrix for the FEM poisson problem should look like at the upper or lower z plane (CWBC = conducting wall boundary condition, TCBC = target corner boundary condition) Screenshot 2024-11-15 at 4 36 52 PM