405: Generic Operator: 1D Nonlinear Poisson equation
Solve the nonlinear Poisson equation
\[-\nabla \varepsilon \nabla u + e^{u}-e^{-u} = f\]
in $\Omega=(0,1)$ with boundary condition $u(0)=0$ and $u(1)=1$ with
\[f(x)= \begin{cases} 1&,x>0.5\\ -1&, x<0.5 \end{cases}.\]
This stationary problem is an example of a nonlinear Poisson equation or Poisson-Boltzmann equation. Such equation occur e.g. in simulations of electrochemical systems and semicondutor devices.
module Example405_GenericOperator
using Printf
using VoronoiFVM
function main(;n=10,Plotter=nothing,verbose=false, unknown_storage=:sparse)
# Create a one-dimensional discretization
h=1.0/convert(Float64,n)
X=collect(0:h:1)
grid=VoronoiFVM.Grid(X)
# A parameter which is "passed" to the flux function via scope
ϵ=1.0e-3
# This generic operator works on the full solution seen as linear vector, and indexing
# into it needs to be performed with the help of idx (defined below for a solution vector)
# Here, instead of the flux function we provide a "generic operator"
# which provides the stiffness part of the problem. Its sparsity is detected automatically
# using SparsityDetection.jl
function generic_operator!(f,u,sys)
f.=0.0
for i=1:length(X)-1
du=ϵ*(u[idx[1,i]]-u[idx[1,i+1]])/(X[i+1]-X[i])
f[idx[1,i]]+=du
f[idx[1,i+1]]-=du
end
end
# Source term
function source!(f,node)
if node[1]<=0.5
f[1]=1
else
f[1]=-1
end
end
# Reaction term
function reaction!(f,u,node)
f[1]=exp(u[1]) - exp(-u[1])
end
# Create a physics structure
physics=VoronoiFVM.Physics(
generic=generic_operator!,
source=source!,
reaction=reaction!)
# Create a finite volume system - either
# in the dense or the sparse version.
# The difference is in the way the solution object
# is stored - as dense or as sparse matrix
sys=VoronoiFVM.System(grid,physics,unknown_storage=unknown_storage)
# Add species 1 to region 1
enable_species!(sys,1,[1])
# Set boundary conditions
boundary_dirichlet!(sys,1,1,0.0)
boundary_dirichlet!(sys,1,2,1.0)
# Create a solution array
inival=unknowns(sys,inival=0.5)
solution=unknowns(sys)
idx=unknown_indices(solution)
# Create solver control info
control=VoronoiFVM.NewtonControl()
control.verbose=verbose
# Stationary solution of the problem
solve!(solution,inival,sys, control=control)
if isplots(Plotter)
coord=coordinates(grid)
Plotter.plot(coord[1,:],solution[1,:],
label="",
title="Nonlinear Poisson",
grid=true,show=true)
end
return sum(solution)
end
function test()
testval=1.5247901344230088
main(unknown_storage=:sparse) ≈ testval && main(unknown_storage=:dense) ≈ testval
end
end
This page was generated using Literate.jl.