Add support for custom boundary conditions

[killah-t-cell]

Now, for a given model, a user can also add custom boundary conditions like so:

Es, xs, vs, t = Plasma.get_vars(plasma; dim=1)
ic = [Es[1](0, xs...) ~ g(xs...)]
bc = [Es[1](t,0.4) ~ f(t)]

sol = Plasma.solve(plasma, dim=1, GPU=false, conditions = [ic;bc]) 

This gives the interface more flexibility to support different simulations.

Additionally, I also added boundary conditions to E in accordance to this

So now, in addition to the boundary and initial conditions every model has to have:

vlasov_ics = [fs[i](time_lb,xs...,vs...) ~ Ps[i](xs,vs) * geometry(xs) for i in eachindex(fs)]
div_E_ic = div_E0 ~ sum([qs[i] * _I(fs[i](time_lb,xs...,vs...)) for i in eachindex(qs)])/ϵ_0 * geometry(xs)
E_bcs = [get_boundaries(t, xs, Es, lb);get_boundaries(t, xs, Es, ub)]

one can also add custom ones.

[killah-t-cell]

cc @SebastianM-C

I think this starts to address what we talked about.

[killah-t-cell]

Refactored a bit. Now we have the ability to add custom boundaries to E, B, and f, but I didn't add this to the interface, as this is more of an advanced use case.

But one could this like so

"""
Get plasma variables to generate custom boundary and initial conditions.

Use get_vars to get variables to set custom initial and boundary conditions like so:

fs, Es, xs, vs, t = Plasma.get_vars(plasma; dim=1)
ic = [Es[1](0, xs...) ~ cos(xs[1])]
bc = [Es[1](t,0.4) ~ cos(t)]

sol = Plasma.solve(plasma, dim=1, GPU=false, E_bcs = [Dirichlet(t, xs, Es, lb);ic;bc]) 

This is only advised for more advanced models with custom currents. For simpler models, set a custom boundary for E with E_bcs or initial condition for f.
"""

function get_vars(plasma; dim=3)
    # get E
    Es = Symbolics.variables(:E, 1:dim; T=SymbolicUtils.FnType{Tuple,Real})

    # get fs
    dis = plasma.distributions
    species = [d.species for d in dis]
    fs = Symbolics.variables(:f, eachindex(species); T=SymbolicUtils.FnType{Tuple,Real})

    # get parameters
    @parameters t
    xs,vs = Symbolics.variables(:x, 1:dim), Symbolics.variables(:v, 1:dim)

    # return values
    if plasma isa ElectrostaticPlasma
        return fs, Es, xs, vs, t
    else 
        Bs = Symbolics.variables(:B, 1:dim; T=SymbolicUtils.FnType{Tuple,Real})
        return fs, Es, Bs, xs, vs, t
    end
end

[SebastianM-C]

Hi, sorry for the delay. I'm looking now over the code and I'm a bit confused by the above E_bcs = [Dirichlet(t, xs, Es, lb);ic;bc] example. Shouldn't it be more like E_bcs = Dirichlet or E_bcs = Neumann (which is the default in the code)?

I tried to test this with

TD = 30000
D = species.D
D_D = Distribution(Maxwellian(TD, D.m), D)
G = Geometry()
plasma = ElectrostaticPlasma([D_D], G)

fs, Es, xs, vs, t = get_vars(plasma; dim = 1)
ic = [Es[1](0, xs...) ~ cos(xs[1])]
bc = [Es[1](t, 0.4) ~ cos(t)]

sol = Plasma.solve(plasma, dim = 1, GPU = false, E_bcs = Dirichlet)

but I'm getting some NaNs showing up in the pde_losses and bcs_losses, but I'm not sure to what does it correspond.

Current loss is: NaN
pde_losses: [3.356915290562938e11, NaN]
bcs_losses: [5.334357060395836, NaN, 0.11948624599487331, 0.14394644995372455, 0.0, 0.0]
Current loss is: NaN
pde_losses: [3.356915290562938e11, NaN]
bcs_losses: [5.334357060395836, NaN, 0.11948624599487331, 0.14394644995372455, 0.0, 0.0]
Current loss is: NaN
pde_losses: [NaN, NaN]
bcs_losses: [NaN, NaN, NaN, NaN, NaN, NaN]

I'm not sure if I'm not missing something on my side, I'll let you know if I find anything.

Also, as a diagnostic, I usually just simulate a laser without the plasma before simulating the whole system in order to check for correctness. Do you think that would be possible as an edge case (i.e. no plasma and just EM wave propagation)? I was wondering if this would be useful to debug the above NaNs (by looking just at the external electric field).

[killah-t-cell]

@SebastianM-C the NaN is due to an error in NeuralPDE (an upstream package of Plasma.jl). I will fix the problem in NeuralPDE and then merge this branch (this is the only thing holding it up).

Shouldn't it be more like E_bcs = Dirichlet or E_bcs = Neumann

It is like that now. I refactored the code after I left the comment.

Also, as a diagnostic, I usually just simulate a laser without the plasma before simulating the whole system in order to check for correctness. Do you think that would be possible as an edge case (i.e. no plasma and just EM wave propagation)? I was wondering if this would be useful to debug the above NaNs (by looking just at the external electric field).

We can support that in the future (feel free to open an issue), but for now I will try to just get the benchmarks working.

[fluffynukeit]

Hi. I am interested in trying out this package, but I am having trouble understanding the code and examples. It looks like the naming convention for the equations and variables are lifted out of some physics text or document, from which one excerpt is provided above. I get E and B, but it is not 100% obvious to me what beta, f, x, and v are, or assumptions about their form or units. Can you document this information in the readme or other appropriate location?

[killah-t-cell]

Hi @fluffynukeit, what particularly are you having trouble understanding? I am happy to answer more specific questions or help you with your code if you share it. Feel free to open an issue. The beta in that case just means a certain species of particle.

And f is the distribution, x is the configuration space coordinates and v is the velocity space coordinates.

[fluffynukeit]

I would appreciate having proper definitions for the terms used in the equations, or a link to the reference material you are sourcing from. For instance, f is the "distribution". Distribution of what? I'm guessing here, so maybe something like "f_beta is the probability density function for a particle of species beta at 3D configuration coordinates x and 3D velocity v." I don't know if that's correct for f_beta or not, but is along the lines of specificity that I think would help in understanding the code.

[killah-t-cell]

Take a look at our read me and at this here https://en.m.wikipedia.org/wiki/Vlasov_equation

it will give you some information about Plasma physics. Make sure to also check out PlasmaBenchmarks.jl (that is still a WIP, but there are a few examples there we will use as benchmarks for the package)