Scilab Website | Contribute with GitLab | Mailing list archives | ATOMS toolboxes
Scilab Online Help
2025.0.0 - English


SUNDIALS Linear Solvers

SUNDIALS direct and iterative, preconditionned linear solvers

Syntax

... = solver( ... , linearSolver = linSol, precType = side, precBand = [mu,ml])

Linear solvers used by SUNDIALS solvers

SUNDIALS solvers may use underlying linear solvers when the selected method needs to solve a linear algebraic system during a method step. This is the case when kinsol is used with the Newton or Picard method or when ODE/DAE solvers use an implicit scheme (arkode uses such a solver with explicit schemes also when a mass matrix is present). When the Jacobian is explicitely provided by the user, either dense, banded or sparse, or implicitely by means of jacBand or jacPattern options, then the corresponding direct solver is used, i.e. Lapack's dense or band solver or KLU sparse solver. The default behavior when a linear solver is needed is to use the dense solver and difference quotients to approximate the Jacobian.

When very large scale problems are considered, matrix-free, scaled, preconditioned, iterative linear solvers are sometimes the only way to go. Although SUNDIALS allows to use this class of solver for ODEs, DAEs and nonlinear equations, for the time being the Scilab module only allows them for ODEs (cvode and arkode solvers) and built-in band preconditioning. The related options are the following:

linearSolver

The linear solver given as a string:

  • "CG" : Conjugate Gradient (system matrix must be symmetric)
  • "BCGS" : Stabilized Bi-Conjugate Gradient
  • "FGMR" : Flexible Generalized Minimum Residual
  • "GMR" : Generalized Minimum Residual
  • "TFQMR" : Transpose-Free Qusi-Minimum Residual

linSolMaxIters

The maximum number of iterations/Krylov vectors

precType The preconditioning type or side: "RIGHT" (the default), "LEFT", "BOTH" or "NONE".
precBand

A vector of two positive integers [mu,ml] giving respectively the upper and lower half bandwidth of the preconditionner.

Example

An ODE system is generated from the following 2-species diurnal kinetics advection-diffusion PDE system in 2 space dimensions

\begin{array}{rcl} \frac{\partial c_i}{\partial t} &=& K_h \frac{\partial^2 c_i}{\partial x^2} + V \frac{\partial c_i}{\partial x} 
                     + \frac{\partial y}{\partial y}\left( K_v(y)\frac{\partial c_i}{\partial y}\right) + R_i(c_1,c_2,t),\quad i=1,2,\\
                     R_1(c_1,c_2,t) &=& -q_1 c_1 c_3 - q_2 c_1 c_2 + 2 q_3(t) c_3 + q_4(t) c_2,\\
                     R_2(c_1,c_2,t) &=&  q_1 c_1 c_3 - q_2 c_1 c_2 - q_4(t) c_2,\\
                     K_v(y) &=& K_{v_0} \exp(y/5),
                 \end{array}

K_h, V, K_{v_0}, q_1, q_2, c_3 are constants and q_3(t) and q_4(t) vary diurnaly. The problem is posed on the square [0,20] x [30, 50] (in km), with homogeneous Neumann boundary conditions, and for time t in [0,129600] (three days and a half). The PDE system is treated by central differences on a uniform 100 x 100 mesh, with simple polynomial initial profiles.

The problem is solved with cvode, with the BDF/GMRES method (i.e. using the SUNLinSol_SPGMR linear solver) by passing option linearSolver="GMR" and a banded preconditioner, generated by difference quotients, using the option precBand. The problem is solved with left and right preconditioning (this is the default).

In the example below the right-hand side is a compiled C-gateway from the examples of SUNDIALS library. The Scilab interface to cvode allows to manipulate the state of the discretized PDE as a (2,100,100) 3-dimensional array and the output C of the solver is a (2,100,100,n) 4-dimensional array where n is the length of the tspan vector. A callback is used to represent the two species with a pseudo-color plot.

function out=cbFn(t, C, flag, stats)
    global h1 h2
    out = %f;
    if flag == "init"
        clf
        gcf().axes_size = [740 300];
        gcf().figure_name = "2-species diurnal kinetics advection-diffusion 2D PDE"
        gcf().color_map = parula(128);
        subplot(1,2,1)
        title("Species 1")
        Sgrayplot(x,y,squeeze(C(1,:,:)),zminmax=[0 51146186])
        h1 = gce().children;
        colorbar;
        subplot(1,2,2)
        Sgrayplot(x,y,squeeze(C(2,:,:)),zminmax=[0 1.074D+12])
        h2 = gce().children;
        colorbar;
        title("Species 2")
    elseif flag == "step"
        drawlater
        h1.data(:,3)=C(1,:)';
        h2.data(:,3)=C(2,:)';
        drawnow
        gcf().info_message=sprintf("time = %g s",t);
    end
end

source = fullfile(SCI,"modules","differential_equations","tests","unit_tests","cvDiurnal_kry_bp.c")
SUN_Clink("cvDiurnal_kry_bp_rhs",source,load=%t);

PI = 3.1415926535898;
TWOHR = 7200.0      // number of seconds in two hours
HALFDAY = 4.32e4    // number of seconds in a half day
KH = 4.0e-6         // horizontal diffusivity Kh
VEL = 0.001         // advection velocity V
KV0 = 1.0e-8        // coefficient in Kv(y)

RTOL = 1e-5;
ATOL = 100*RTOL;

MX = 100; // MX and MY values are hard coded in cvDiurnal_kry_bp.c
MY = 100; // don't change these values !
XMIN = 0; XMAX = 20;
YMIN = 30; YMAX = 50;
DX = (XMAX-XMIN)/(MX-1);
DY = (YMAX-YMIN)/(MY-1);
x = XMIN+(0:MX-1)*DX;
y = YMIN+(0:MY-1)*DY;
[X,Y] = ndgrid(x,y);
q4 = 0;
om = PI/HALFDAY;
hdco = KH/DX^2;
haco = VEL/(2*DX);
vdco = KV0/DY^2;
user_data = [q4,om,DX,DY,hdco,haco,vdco];

// initial profiles
C0 = zeros(2,MX,MY);
alph = (0.1*X-1).^2;
alph = 1-alph+0.5*alph.^2;
bet = (0.1*Y-4).^2;
bet = 1-bet+0.5*bet.^2;
C0(1,:,:) = 1e6*alph.*bet;
C0(2,:,:) = 1e12*alph.*bet;
tspan = 0:TWOHR/60:3*HALFDAY;

[t,C,info] = cvode(list("cvDiurnal_kry_bp_rhs",user_data), tspan, C0, t0=0, rtol=RTOL, atol=ATOL,...
    method="BDF", linearSolver="GMR", precBand=[2 2], callback=cbFn);

See also

  • arkode — SUNDIALS ordinary differential equation additive Runge-Kutta solver
  • cvode — SUNDIALS ordinary differential equation solver
  • ida — SUNDIALS differential-algebraic equation solver
  • kinsol — SUNDIALS differential-algebraic equation solver
Report an issue
<< Options (kinsol) Options, features and user functions Options (ODE and DAE solvers) >>

Copyright (c) 2022-2024 (Dassault Systèmes)
Copyright (c) 2017-2022 (ESI Group)
Copyright (c) 2011-2017 (Scilab Enterprises)
Copyright (c) 1989-2012 (INRIA)
Copyright (c) 1989-2007 (ENPC)
with contributors
Last updated:
Thu Oct 24 11:13:06 CEST 2024