Public and private APIs

pdepe(m, pdefunction, icfunction, bdfunction, xmesh, tspan ; params=SkeelBerzins.Params(), kwargs...)

Solve 1D elliptic and parabolic partial differential equation(s) using the spatial discretization method described in [1]. Note that to use this method, one of the PDEs must be parabolic. The time discretization is either done by the implicit Euler method (internal method) or by using a ODE/DAE solver from the DifferentialEquations.jl package. For more information on how to define the different inputs to solve a problem, look at the following sections: Problem Definition and Solvers.

Input arguments:

• m: scalar referring to the symmetry of the problem. It can either take the value m=0, m=1 or m=2 representing cartesian, cylindrical or spherical coordinates respectively.
• pdefunction: Function. Defines the formulation of the PDE(s) using capacity, flux and source terms.
• icfunction: Function. Defines the initial condition(s) of the problem to solve (if tstep!=Inf - initial condition(s) from the ODE/DAE problem, else if tstep=Inf - initial value(s) used for the newton solver).
• bdfunction: Function. Defines the boundary conditions of the problem.
• xmesh: 1D array representing the spatial mesh on which the user intends to obtain the solution.
• tspan: tuple $(t_0, t_{end})$ representing the time interval of the problem.

Keyword argument:

• params: defines a SkeelBerzins.Params struct containing the keyword arguments from the solvers.
• kwargs: instead of using the SkeelBerzins.Params struct, the user can pass directly fields from this particular struct to the solver.

Returns a RecursiveArrayTools.DiffEqArray or a SkeelBerzins.ProblemDefinition struct, depending on the selected solver (either :euler or :DiffEq). The obtained solution includes a linear interpolation method in time and can be used to evaluate the solution at any time step within the interval $(t_0,t_{end})$ (accessible using sol(t)). A spatial interpolation similar as the pdeval function is available on the solution object using the command sol(x_eval,t,pb).

pdepe(m, pdefunction, icfunction, bdfunction, xmesh ; params=SkeelBerzins.Params(tstep=Inf), kwargs...)

Solve 1D elliptic PDE(s) using the spatial discretization method described in [1] - pdepe variant to solve stationary problems. Performs one step of the implicit Euler method.

Input arguments:

• m: scalar referring to the symmetry of the problem. It can either take the value m=0, m=1 or m=2 representing cartesian, cylindrical or spherical coordinates respectively.
• pdefunction: Function. Defines the formulation of the PDE(s) using capacity, flux and source terms (capacity term should be set to 0).
• icfunction: Function. It defines the initial value(s) used for the Newton solver.
• bdfunction: Function. Defines the boundary conditions of the problem.
• xmesh: 1D array representing the spatial mesh on which the user intends to obtain the solution.

Keyword argument:

• params: defines a SkeelBerzins.Params structure containing the keyword arguments from the solvers.
• kwargs: instead of using the SkeelBerzins.Params struct, the user can pass directly fields from this particular struct to the solver.

Returns a 1D Array with the solution available at the points defined by the spatial discretization xmesh.

struct Params

Structure containing all the keyword arguments for the solver pdepe.

• solver::Symbol: Choice of the time discretization either use :euler for internal implicit Euler method or :DiffEq for the DifferentialEquations.jl package.

• tstep::Union{Float64, Vector{Float64}}: Defines a time step (either pass a Float64 or a Vector) when using the implicit Euler method. When set to tstep=Inf, it solves the stationary version of the problem.

• hist::Bool: Flag, returns with the solution, a list of 1d-array with the history from the newton solver.

• sparsity::Symbol: Choice of the type of matrix (:sparseArrays, :banded) use to store the jacobian.

• linsolve::Union{Nothing, LinearSolve.SciMLLinearSolveAlgorithm}: Choice of the solver for the LSE in the newton method, see LinearSolve.jl.

• maxit::Int64: Maximum number of iterations for the Newton solver.

• tol::Float64: Tolerance used for Newton method. Returns solution if $||\; u_{i+1} - u_{i} \;||_2 <$ tol.

• data::Bool: Returns the data of the PDE problem

• nb_design_var::Int64: Number of design variables for PDE Constrained Optimization

ODEProblem(problem)

Generate an ODEProblem from the ODEFunction which then can be solved by using the solve method.

Input arguments:

• problem: Structure of type SkeelBerzins.ProblemDefinition.

reshape(sol_ODE, problem)

Reshape the solution sol_ODE obtained by the DifferentialEquations.jl package to obtain the solution at time t as a 2D-Array of size (npde,Nx).

Indeed since in the spatial discretization, we flattened the problem, the solution has a similar size. So by reshaping the solution, we get a solution organised by unknows.

pdeval(m, xmesh, u_approx, x_eval, pb)

Function that interpolates with respect to the space component the solution obtained by the solver function pdepe.

Input arguments:

• m: symmetry of the problem (m=0,1,2 for cartesian, cylindrical or spherical).
• xmesh: space discretization.
• u_approx: approximate solution obtained by the solver pdepe.
• x_eval: point or vector of points where to interpolate the approximate solution.
• pb: structure defining the problem definition, see SkeelBerzins.ProblemDefinition.

Returns a tuple $(u,dudx)$ corresponding to the solution and its partial derivative with respect to the space component evaluated in x_eval.

assemble!(du, u, problem, t)

Performs space discretization following the difference equations described in [1].

Assemble the right-hand side $f$ to generate an ODE/DAE problem:

\[(M) \frac{du}{dt} = f(u,problem,t)\]

where the input problem is defined as a SkeelBerzins.ProblemDefinition structure.

This function is specified in a way that it is compatible with the DifferentialEquations.jl package.

interpolation(xl, ul, xr, ur, quadrature_point, problem)

Interpolate u and $\frac{du}{dx}$ between two discretization points at some specific quadrature point. Method use for scalar PDE.

Input arguments:

• xl: left boundary of the current interval.
• ul: solution evaluated at the left boundary of the current interval.
• xr: right boundary of the current interval.
• ur: solution evaluated at the right boundary of the current interval.
• quadrature_point: quadrature point chosen according to the method described in [1].
• problem: Structure of type SkeelBerzins.ProblemDefinition.

interpolation(xl, ul, xr, ur, quadrature_point, m, singular, npde)

Interpolate u and $\frac{du}{dx}$ between two discretization points at some specific quadrature point and return them as static vectors. Method use for system of PDEs.

Input arguments:

• xl: left boundary of the current interval.
• ul: solution evaluated at the left boundary of the current interval.
• xr: right boundary of the current interval.
• ur: solution evaluated at the right boundary of the current interval.
• quadrature_point: quadrature point chosen according to the method described in [1].
• m: symmetry of the problem (given as a type).
• singular: indicates whether the problem is regular or singular (given as a type).
• npde: number of PDEs (given as a type).

newton(b, tau, timeStep, problem, mass_matrix, cache, rhs ; tol=1.0e-10, maxit=100, hist_flag=false, linSol=nothing)

Newton method solving nonlinear system of equations. The Jacobi matrix used for the iteration rule is computed with the help of the SparseDiffTools.jl package.

Input arguments:

• b: right-hand side of the system to solve.
• tau: constant time step used for the time discretization.
• timeStep: current time step of tspan.
• problem: Structure of type SkeelBerzins.ProblemDefinition.
• mass_matrix: mass matrix of the problem, see [SkeelBerzins.mass_matrix][@ref].
• cache: SparseDiffTools.ForwardColorCache. To avoid allocating the cache in each iteration of the newton solver when computing the jacobian.
• rhs: preallocated vector to avoid creating allocations.

Keyword arguments:

• tol: tolerance or stoppping criteria (by default to 1.0e-10).
• maxit: maximum number of iterations (by default to 100).
• hist_flag: flag to save the history and returns it (by default to false).
• linSol: choice of the solver for the LSE, see LinearSolve.jl (by default nothing).

Returns the solution of the nonlinear system of equations and if hist_flag=true, the history of the solver.

newton_stat(b, tau, timeStep, problem, cache, rhs ; tol=1.0e-10, maxit=100, hist_flag=false, linSol=nothing)

Newton method solving nonlinear system of equations (variant of newton for stationary problems).

struct ProblemDefinition{T1, T2, T3, Tv<:(AbstractVector), Ti<:Integer, Tm<:Number, elTv<:Number, pdeFunction<:Function, icFunction<:Function, bdFunction<:Function}

Structure storing the problem definition.

• npde::Integer: Number of unknowns

• Nx::Integer: Number of discretization points

• xmesh::AbstractVector: Grid of the problem

• tspan::Tuple{Tm, Tm} where Tm<:Number: Time interval

• singular::Bool: Flag to know if the problem is singular or not

• m::Integer: Symmetry of the problem

• jac::Union{BandedMatrices.BandedMatrix{elTv, Matrix{elTv}, Base.OneTo{Ti}}, SparseArrays.SparseMatrixCSC{elTv, Ti}} where {Ti<:Integer, elTv<:Number}: Jacobi matrix

• nb_design_var::Integer: Number of design variables for PDE Constrained Optimization

• inival::Vector{elTv} where elTv<:Number: Evaluation of the initial condition

• ξ::Vector{elTv} where elTv<:Number: Interpolation points from the paper

• ζ::Vector{elTv} where elTv<:Number

• pdefunction::Function: Function defining the coefficients of the PDE

• icfunction::Function: Function defining the initial condition

• bdfunction::Function: Function defining the boundary conditions

problem_init(m, xmesh, tspan, pdefun, icfun, bdfun, params)

Function initializing the problem.

Input arguments: similar as pdepe.

Returns the size of the space discretization Nx, the number of PDEs npde, the initial value inival, some data types elTv and Ti, and the struct containing the problem definition pb.

get_sparsity_pattern(sparsity, Nx, npde, elTv)

Function that provides the sparsity pattern in a SparseMatrixCSC.

get_sparsity_pattern(sparsity, Nx, npde, elTv)

Function that provides the sparsity pattern in a BandedMatrix.

get_quad_points_weights(m, alpha, beta, gamma, singular)

Calculate the quadrature points and weights for the one-point Gauss quadrature based on the problem's specific symmetry, as described in the paper [1].

Input arguments:

• m: scalar representing the symmetry of the problem.
• alpha: 1D array containing the left boundaries of the subintervals.
• beta: 1D array containing the right boundaries of the subintervals.
• gamma: 1D array containing the middle points of the subintervals.
• singular: boolean indicating whether the problem is singular or not.

Returns the quadrature points xi and the weights zeta.

implicitEuler!(y,u,problem,tau,mass_matrix,timeStep)

Assemble the system for the implicit Euler method.

implicitEuler_stat!(y,u,problem,tau,timeStep)

Assemble the system for the implicit Euler method (variant method for stationary problems).

mass_matrix(problem)

Assemble the diagonal mass matrix M of the system of differential equations when solving a problem with at least one parabolic PDE. The coefficients from M either take the value 0 or 1 since it is scaled in the difference equations in the right-hand side.

The entries of the matrix are set to 0 when the corresponding equation of the system is elliptic or the boundary condition is pure Dirichlet leading to solve a Differential-Algebraic system of Equations. In the case where the mass matrix is identity, we solve a system of ODEs.

interpolate_sol_time(u_approx,t)

Linear interpolatation of the solution with respect to the time component.

Input arguments:

• u_approx: approximate solution obtained by the solver pdepe.
• t: time $t \in [t_0,t_{end}]$.

interpolate_sol_space(u_approx, x_eval, pb, times ; val=true, deriv=true)

Similar as pdeval, i.e. interpolate the solution and its derivative with respect to the space component.

Input arguments:

• u_approx: approximate solution obtained by the solver pdepe.
• x_eval: point or vector of points where to interpolate the approximate solution.
• times: point or vector of points which denotes the time step of the solution.
• pb: structure defining the problem definition, see SkeelBerzins.ProblemDefinition.

Keyword arguments:

• val: flag to return the evaluation of the solution at the x_eval position.
• deriv: flag to return the derivative with respect to the space component of the solution at the x_eval position.

Returns a tuple $(u,dudx)$ corresponding to the solution and its partial derivative with respect to the space component evaluated in x_eval at time times.

ODEFunction(problem)

Generate an ODEFunction from the spatial discretization derived in the SkeelBerzins.assemble! function. It is expressed as a mass matrix ODE if the mass matrix is different from the identity matrix or as a simple system of ODEs otherwise and defines the problem with respect to the sparsity pattern.

Input argument:

• problem: Structure of type SkeelBerzins.ProblemDefinition.