Coupling with Nernst-Planck-Poisson

[j-fu]

As discussed, I try the coupling with Nernst-Planck-Poisson.
First thing to try is a pore with inlet and outlet

                
                wall
     +-----------------------+
in   |                       | out
     |                       | 
     +-----------------------+
                wall

So I have some questions:

• How to set "do nothing boundary conditions" at inlet and outlet ?
• In the C++ pdelib, we had Femsolver.AssignElectricFieldAndCharges(potential, charge_density). Is is possible to have this here as well ? Or is there a way to specify an element-wise body force (which I could calculate if necessary for the time being)?

[chmerdon]

For the "do nothing boundary conditions" you have to do nothing: if no operators for boundary conditions are specified all degrees of freedom on the boundary are kept free and if the Neumann data is zero, nothing needs to be done. I guess 'wall' has an dedicated boundary region number r. Then, you need only assign boundary conditions there:

bnd_op_wall = HomogeneousBoundaryData(u; regions = [r])

which you add to the ProblemDescription as usual via

assign_operator!(PD, bnd_op_wall)

---

Concerning the Lorentz force:
In case the grids match in the sense that the FVM solution gives values for the nodes, I would suggest to first construct an FEVector that can carry everything: velocity, pressure, electrostatic potential, concentrations. For this we need finite element spaces:

FES_u = FESpace{H1BR{2}}(grid) # Bernardi-Raugel velocity space
FES_p = FESpace{L2P0{1}}(grid) # P0 pressure space appropriate for BR above
FES_V = FESpace{H1P1{1}}(grid)
FES_c = FESpace{H1P1{nspecies}}(grid)

sol = FEVector([FES_u, FES_p, FES_V, FES_c])

Then we can write a kernel for the Lorentz force, probably something like

function lorentz_kernel(result, args, qpinfo)
params = qpinfo.params
q = sum(params .* view(args, 3:nspecies+2)) # = total charge
result .= q * view(args,1:2) # = q \nabla V
end

The operator would look like

LinearForm(lorentz_kernel, [id(1)], [grad(3), id(4)]; params = charge_numbers)

Therein, the numbers correspond to the number of the FESpaces above, so grad(3) gives the gradient of the electrostatic potential and id(4) the concentrations of the species, which enter the kernel in that order via the argument args above.
Via the params argument constant coefficients can be funneled into the kernel, e.g. the charge numbers of the species in this case.

This LinearForm is then added to the ProblemDescription of a classical Stokes problem for the unknowns u and p.

[chmerdon]

ah, and in each iteration the coefficients for V and c can be overwritten by

sol[3] .= new_coefficients_for_V
sol[4] .= new_coefficients_for_c

Update: If you also create Unknowns [u,p,V,c] (via u = Unknown("velocity") etc.) and assign them in the FEVector constructor via (; tags = [u,p,V,c]), you can use u,p,V,c instead of 1,2,3,4 everywhere.

[chmerdon]

And the solve call for the Stokes problem (assumed to have the ProblemDescription PD) looks a bit different:

sol = solve(PD, [FES_u, FES_p]; init = sol)

such that only the u and p blocks are solved (but the others can be used for assembly). I hope this works, we can fiddle out the details on Monday, if you like.