Skip to content

Un solveur multigrille géométrique pour Mathematica ?

Si vous rencontrez une partie que vous ne comprenez pas, vous pouvez nous laisser un commentaire et nous vous répondrons dans les plus brefs délais.

Solution :

Fond d'écran :

Des détails sur les solveurs multigrilles peuvent être trouvés dans ce script assez soigné de Volker John. C'est essentiellement la source dans laquelle j'ai puisé les informations pour mettre en œuvre le solveur de cycle en V ci-dessous.

En un mot, un solveur multigrille s'appuie sur deux ingrédients : Une hiérarchie de systèmes linéaires (avec ce qu'on appelle opérateurs de prolongation entre eux) et une famille de lisseurs.

Implémentation

Méthode CG comme lisseur

J'utilise un solveur itératif de gradient conjugué comme lisseur. Pour une raison quelconque, le solveur de gradient conjugué de Mathematica a une latence exorbitante, donc j'utilise ma propre implémentation que j'ai écrite il y a plusieurs années. Elle est très facile à mettre en œuvre ; tous les détails nécessaires peuvent être trouvés, par exemple, ici. Notez que mon implémentation renvoie un Association qui fournit également des informations sur le processus de résolution. (En particulier pour les EDP transitoires à coefficients variables, le nombre d'itérations nécessaires pour réduire le résidu en dessous d'une tolérance donnée est souvent une information précieuse que l'on peut vouloir utiliser, par exemple, pour déterminer si le préconditionneur doit être mis à jour).

ClearAll[CGLinearSolve];

Options[CGLinearSolve] = {
   "Tolerance" -> 10^(-8),
   "StartingVector" -> Automatic,
   MaxIterations -> 1000,
   "Preconditioner" -> Identity
   };

CGLinearSolve[A_?SquareMatrixQ, b_?MatrixQ, opts : OptionsPattern[]] :=
   CGLinearSolve[A, #, opts] & /@ Transpose[b]

CGLinearSolve[A_?SquareMatrixQ, b_?VectorQ, OptionsPattern[]] :=

  Module[{r, u, δ, δ0, p, ρ, ρold, 
    z, α, β, x, TOL, iter, P, precdata, normb, maxiter},
   P = OptionValue["Preconditioner"];
   maxiter = OptionValue[MaxIterations];
   normb = Sqrt[b.b];
   If[Head[P] === String,
    precdata = SparseArray`SparseMatrixILU[A, "Method" -> P];
    P = x [Function] SparseArray`SparseMatrixApplyILU[precdata, x];
    ];
   If[P === Automatic,
    precdata = SparseArray`SparseMatrixILU[A, "Method" -> "ILU0"];
    P = x [Function] SparseArray`SparseMatrixApplyILU[precdata, x];
    ];
   TOL = normb OptionValue["Tolerance"];
   If[OptionValue["StartingVector"] === Automatic,
    x = ConstantArray[0., Dimensions[A][[2]]];
    r = b
    ,
    x = OptionValue["StartingVector"];
    r = b - A.x;
    ];
   z = P[r];
   p = z;
   ρ = r.z;
   δ0 = δ = Sqrt[r.r];
   iter = 0;
   While[δ > TOL && iter < maxiter,
    iter++;
    u = A.p;
    α = ρ/(p.u);
    x = x + α p;
    ρold = ρ;
    r = r - α u;
    δ = Sqrt[r.r];
    z = P[r];
    ρ = r.z;
    β = ρ/ρold;
    p = z + β p;
    ];
   Association[
    "Solution" -> x,
    "Iterations" -> iter,
    "Residual" -> δ,
    "RelativeResidual" -> Quiet[Check[δ/δ0, ∞]],
    "NormalizedResidual" -> Quiet[Check[δ/normb, ∞]]
    ]
   ];

Lissage de Jacobi pondéré

La méthode de Jacobi pondérée est un solveur itératif très simple qui emploie des itérations de Richardson avec la diagonale de la matrice comme préconditionneur (la matrice ne doit pas avoir d'éléments nuls sur la diagonale !) et un peu d'amortissement. Ne fonctionne en général que pour les matrices à dominante diagonale et les matrices à définition positive (si la matrice "Weight" est choisi suffisamment petit). Ce n'est pas vraiment mauvais mais ce n'est pas excellent non plus. Dans le problème test ci-dessous, elle nécessite généralement un cycle V de plus comme lisseur CG.

Options[JacobiLinearSolve] = {
   "Tolerance" -> 10^(-8),
   "StartingVector" -> Automatic,
   MaxIterations -> 1000,
   "Weight" -> 2./3.
   };

JacobiLinearSolve[A_?SquareMatrixQ, b_?MatrixQ, opts : OptionsPattern[]] := 
  JacobiLinearSolve[A, #, opts] & /@ Transpose[b]

JacobiLinearSolve[A_?SquareMatrixQ, b_?VectorQ, OptionsPattern[]] := 
  Module[{ω, x, r, ωd, dd, iter, δ, δ0, normb, TOL, maxiter},
   ω = OptionValue["Weight"];
   maxiter = OptionValue[MaxIterations];
   normb = Max[Abs[b]];
   TOL = normb OptionValue["Tolerance"];
   If[OptionValue["StartingVector"] === Automatic,
    x = ConstantArray[0., Dimensions[A][[2]]];
    r = b;
    ,
    x = OptionValue["StartingVector"];
    r = b - A.x;
    ];
   ωd = ω/Normal[Diagonal[A]];
   δ = δ0 = Max[Abs[r]];
   iter = 0;
   While[δ > TOL && iter < maxiter,
    iter++;
    x += ωd r;
    r = (b - A.x);
    δ = Max[Abs[r]];
    ];
   Association[
    "Solution" -> x,
    "Iterations" -> iter,
    "Residual" -> δ,
    "RelativeResidual" -> Quiet[Check[δ/δ0, ∞]],
    "NormalizedResidual" -> Quiet[Check[δ/normb, ∞]]
    ]
   ];

Mise en place du solveur

Ensuite, il y a une fonction qui prend la matrice du système Matrix et une famille d'opérateurs de prologation Prolongations et crée une GMGLinearSolveFunction . Cet objet contient une méthode de résolution linéaire pour le niveau le plus profond de la hiérarchie, les opérateurs de prolongation et - dérivés de Matrix et Prolongations - une matrice de système linéaire pour chaque niveau de la hiérarchie.

Comme c'est l'idée centrale des schémas de Galerkin en FEM, nous interprétons la matrice système Matrix sur la grille la plus fine comme un opérateur linéaire $A_0 colon X_0 to X_0'$$X_0$ désigne l'espace de fonctions d'éléments finis de fonctions continues, linéaires par morceaux, sur la maille la plus fine et $X_0'$ désigne son espace dual. En dénotant l'espace de fonctions d'éléments finis sur la maille $X_0'$. $i$ -ième sous-grille par $X_i$ et en interprétant les opérateurs de prolongation dans la liste Prolongations comme des enchâssements linéaires $J_i colon X_{i} hookrightarrow X_{i-1}$ , on obtient les opérateurs linéaires $A_i colon X_i to X_i'$ par Sous-espace de Galerkin projection pullback, c'est-à-dire en exigeant que le diagramme suivant soit commutatif :

$$require{AMScd}
begin{CD}
X_0 @<{J_1}< X_1 @<{J_2}< dotsm @<{J_{n-1}}< X_{n-1} @<{J_{n}}< X_n\ @VV{A_0}V @VV{A_1}V @. @VV{A_{n-1}}V @VV{A_{n}}V\ X_0' @>{J_1'}> X_1' @>{J_2'}> dotsm @>{J_{n-1}'}> X_{n-1}' @>{J_{n}'}> X_n'
end{CD}$$

Par défaut, LinearSolve est utilisée comme solveur pour la grille la plus grossière, mais l'utilisateur peut spécifier toute autre fonction... F via l'option "CoarsestGridSolver" -> F.

Quelques jolies impressions pour le créé GMGLinearSolveFunction créés est également ajoutée.

ClearAll[GMGLinearSolve, GMGLinearSolveFunction];

GMGLinearSolve[
   Matrix_?SquareMatrixQ,
   Prolongations_List,
   OptionsPattern[{
     "CoarsestGridSolver" -> LinearSolve
     }]
   ] := Module[{A},
   (*Galerkin subspace projections of the system matrix*)
   A = FoldList[Transpose[#2].(#1.#2) &, Matrix, Prolongations];
   GMGLinearSolveFunction[
    Association[
     "MatrixHierarchy" -> A,
     "Prolongations" -> Prolongations,
     "CoarsestGridSolver" -> OptionValue["CoarsestGridSolver"][A[[-1]]],
     "CoarsestGridSolverFunction" -> OptionValue["CoarsestGridSolver"]
     ]
    ]
   ];

GMGLinearSolveFunction /: 
  MakeBoxes[S_GMGLinearSolveFunction, StandardForm] := 
  BoxForm`ArrangeSummaryBox[GMGLinearSolveFunction, "", 
   BoxForm`GenericIcon[LinearSolveFunction],
   {
    {
     BoxForm`MakeSummaryItem[{"Specified elements: ", 
       Length[S[[1, "MatrixHierarchy", 1]]["NonzeroValues"]]}, 
      StandardForm]
     },
    {
     BoxForm`MakeSummaryItem[{"Dimensions: ", 
       Dimensions[S[[1, "MatrixHierarchy", 1]]]}, StandardForm],
     BoxForm`MakeSummaryItem[{"Depth: ", 
       Length[S[[1, "MatrixHierarchy"]]]}, StandardForm]
     }
    },
   {
    BoxForm`MakeSummaryItem[{
      Invisible["Dimensions: "],
      Column[Dimensions /@ S[[1, "MatrixHierarchy", 2 ;;]]]},
     StandardForm],
    BoxForm`MakeSummaryItem[{
      "CoarsestGridSolver: ",
      S[[1, "CoarsestGridSolverFunction"]]
      }, StandardForm]
    },
   StandardForm, "Interpretable" -> False
   ];

Le solveur

Ce qui suit est le solveur réel du cycle en V. À ma propre surprise, l'algorithme n'était pas si difficile à mettre en œuvre. Comme toujours, la plus grande partie du travail a dû être investie dans l'interface utilisateur (et elle n'est toujours pas complète car il manque une protection contre les erreurs 1D-10T).

En fait, ce solveur de cycle en V est un solveur purement algébrique (AMG) ; le géométrie dans "geometric multigrid solver" réside dans la façon dont la hiérarchie matricielle a été construite (à savoir par des grilles géométriquement imbriquées et des méthodes de sous-espace Galerkin).

Options[GMGLinearSolveFunction] = {
   "StartingVector" -> Automatic,
   "Tolerance" -> 1. 10^-8,
   "MaxIterations" -> 25,
   "StartingVectorSmoothingCounts" -> 12,
   "PreSmoothingCounts" -> 8,
   "PostSmoothingCounts" -> 8,
   "Smoother" -> Function[
     {A, b, x0, ν, p},
     (
      CGLinearSolve[A, b,
        MaxIterations -> ν,
        "StartingVector" -> x0,
        "Tolerance" -> 10^-12
        ]["Solution"]
      )
     ],
   "SmootherParameters" -> None
   };

GMGLinearSolveFunction /: GMGLinearSolveFunction[a_Association][
   Rhs_?VectorQ,
   opts0 : OptionsPattern[]
   ] := With[{
    J = a["Prolongations"],
    A = a["MatrixHierarchy"],
    Asol = a["CoarsestGridSolver"]
    },
   Module[{smoother, Rhsnorm, p, n, v, f, depth, allocationtime, startingvector, startingvectortime, solvetime, startingvectorresidual, residual, ν0, ν1, ν2, tol, maxiter, iter, opts},
    opts = Merge[{
       Options[GMGLinearSolveFunction],
       opts0
       }, Last
      ];
    n = Length /@ A;
    depth = Length[n];

    smoother = opts["Smoother"];
    p = opts["SmootherParameters"];
    If[p === None, p = ConstantArray[{}, depth];];

    (* allocate memory for computations *)

    allocationtime = AbsoluteTiming[
       v = ConstantArray[0., #] & /@ n;
       f = Join[{Rhs}, ConstantArray[0., #] & /@ Most[n]];
       ][[1]];

    (* compute starting vector *)

    startingvectortime = AbsoluteTiming[
       If[VectorQ[opts["StartingVector"]],
        v[[1]] = opts["StartingVector"];
        ,
        If[opts["StartingVector"] =!= "Null", opts["StartingVector"] == Automatic];];

       If[opts["StartingVector"] === Automatic,
        Module[{b},
          ν0 = opts["StartingVectorSmoothingCounts"];
          If[! ListQ[ν0], ν0 = If[IntegerQ[ν0], ConstantArray[ν0, Length[n] - 1], ν0 /@ Range[depth]]];
          b = FoldList[#1.#2 &, Rhs, J];
          v[[depth]] = Asol[b[[depth]]];

          Do[v[[i]] = smoother[A[[i]], b[[i]], J[[i]].v[[i + 1]], ν0[[i]], p[[i]]], {i, depth - 1, 1, -1}];
          ];
        ,
        ν0 = None;
        ];
       ][[1]];
    startingvector = v[[1]];
    residual = startingvectorresidual = Max[Abs[Rhs - A[[1]].startingvector]];

    (* perform V-cycles until tolerance is met *)

    solvetime = AbsoluteTiming[
       ν1 = opts["PreSmoothingCounts"];
       If[! ListQ[ν1], ν1 = If[IntegerQ[ν1], ConstantArray[ν1, Length[n] - 1], ν1 /@ Range[depth]]];
       ν2 = opts["PostSmoothingCounts"];
       If[! ListQ[ν2], ν2 = If[IntegerQ[ν2], ConstantArray[ν2, Length[n] - 1], ν2 /@ Range[depth]]];
       Rhsnorm = Max[Abs[Rhs]];
       tol = opts["Tolerance"] Rhsnorm;
       maxiter = opts["MaxIterations"];
       iter = 0;
       While[
        residual > tol && iter < maxiter,
        iter++;
        Do[
         v[[i]] = smoother[A[[i]], f[[i]], N[Boole[i == 1]] v[[i]], ν1[[i]], p[[i]]];
         f[[i + 1]] = (f[[i]] - A[[i]].v[[i]]).J[[i]],
         {i, 1, depth - 1}];

        (* solve at deepest level with "CoarsestGridSolver" *)

         v[[depth]] = Asol[f[[depth]]];

        Do[
         v[[i]] = smoother[A[[i]], f[[i]], v[[i]] + J[[i]].v[[i + 1]], ν2[[i]], p[[i]]],
         {i, depth - 1, 1, -1}];
        residual = Max[Abs[Subtract[Rhs, A[[1]].v[[1]]]]];
        ];
       ][[1]];

    Association[
     "Solution" -> v[[1]],
     "StartingVectorResidual" -> startingvectorresidual,
     "StartingVectorNormalizedResidual" -> 
      startingvectorresidual/Rhsnorm,
     "Residual" -> residual,
     "NormalizedResidual" -> residual/Rhsnorm,
     "SuccessQ" -> residual < tol,
     "Timings" -> [email protected][
        "Total" -> allocationtime + startingvectortime + solvetime,
        "Allocation" -> allocationtime,
        "StartingVector" -> startingvectortime,
        "V-Cycle" -> solvetime
        ],
     "V-CycleCount" -> iter,
     "SmootingCounts" -> [email protected][
        "StartingVector" -> {ν0},
        "Pre" -> {ν1},
        "Post" -> {ν2}
        ],
     "StartingVector" -> startingvector,
     "Smoother" -> smoother,
     "Depth" -> depth
     ]
    ]
   ];

Application

Ce dont nous avons besoin maintenant, c'est d'un cas de test ! Par hasard, j'ai récemment mis à jour ma routine de subdivision de boucle de telle sorte qu'elle renvoie également la matrice de subdivision si nous le demandons gentiment. Nous pouvons utiliser ces matrices de subdivision comme opérateurs de prolongation !

Donc, commençons par un maillage assez grossier sur le disque unitaire et raffinons-le par Loop subdivision (vous aurez besoin du code de LoopSubdivide si vous voulez essayer cela) :

R = DiscretizeRegion[Disk[], MaxCellMeasure -> 0.001];
depth = 5;
{R, J} = {Last[#1], Reverse[Rest[#2]]} & @@ 
   [email protected][LoopSubdivide, {R, {{0}}}, depth - 1];

Résolvons le problème elliptique suivant avec des conditions aux limites de Neumann sur le disque.$varOmega$:

$$begin{array}{rcll}
(varepsilon - Delta) , u &= &f, & text{in $varOmegasetminus partial varOmega$,}\N
nu ,u&= &0, & text{on $varOmega$partiel varOmega$,}
end{array}$$

Nous pouvons utiliser les capacités FEM de Mathematica pour assembler la matrice du système et le côté droit pour nous :

f = X [Function] 
   16. Sinc[4. Pi Sqrt[Abs[Dot[X + 0.5, X + 0.5]]]] - 
    16. Sinc[4. Pi Sqrt[Abs[Dot[X - 0.5, X - 0.5]]]] + 
    N[Sign[X[[2]]]] + N[Sign[X[[1]]]];
fvec = Map[f, MeshCoordinates[R]];

Needs["NDSolve`FEM`"]
Rdiscr = ToElementMesh[
   R,
   "MeshOrder" -> 1,
   "NodeReordering" -> False,
   MeshQualityGoal -> 0
   ];
vd = NDSolve`VariableData[{"DependentVariables", "Space"} -> {{u}, {x, y}}];
sd = NDSolve`SolutionData[{"Space"} -> {Rdiscr}];
cdata = InitializePDECoefficients[vd, sd,
   "DiffusionCoefficients" -> {{-IdentityMatrix[2]}}, 
   "MassCoefficients" -> {{1}}, "LoadCoefficients" -> {{f[{x, y}]}}];
bcdata = InitializeBoundaryConditions[vd, 
   sd, {DirichletCondition[u[x, y] == 0., True]}];
mdata = InitializePDEMethodData[vd, sd];
dpde = DiscretizePDE[cdata, mdata, sd];

{b, L, damping, M} = dpde["All"];
b = Flatten[Normal[b]];
A = L + 0.0001 M;

Maintenant, nous créons une GMGLinearSolveFunction et résolvons l'équation :

S = GMGLinearSolve[A, J]; // AbsoluteTiming // First
xGMG = S[b,
      "Tolerance" -> 1. 10^-8,
      "StartingVectorSmoothingCounts" -> 12,
      "PreSmoothingCounts" -> 8,
      "PostSmoothingCounts" -> 8
      ]["Solution"]; // AbsoluteTiming // First

0.835408

1.04969

Timings

Voici les timings de quelques autres solveurs :

xKrylov = LinearSolve[A, b, Method -> {
       "Krylov",
       "Method" -> "ConjugateGradient",
       "Preconditioner" -> "ILU0"
       }]; // AbsoluteTiming // First
xTAUCS = LinearSolve[A, b, "Method" -> "Cholesky"]; // AbsoluteTiming // First
xUMFPACK = LinearSolve[A, b]; // AbsoluteTiming // First
xPardiso = LinearSolve[A, b, "Method" -> "Pardiso"]; // AbsoluteTiming // First

67.948

6.89134

6.0961

2.30715

Trois choses à observer ici :

  1. Mathematica "ConjugateGradient" est le absolue perdant ici. (Mais ne me demandez pas pour "GMRES" ou "BiCGSTAB"; je n'étais pas d'humeur à les attendre).

  2. "Cholesky" ne peut pas convertir sa limitation aux matrices définies positives en un quelconque avantage. C'est aussi pourquoi je ne l'utilise jamais.

  3. GMGLinearSolve est en fait un peu plus rapide que "Pardiso".

Erreurs

Voici les erreurs ; j'utilise la solution de UMFPACK 's comme "vérité de base" (cela n'a pas d'importance, cependant).

Max[Abs[xUMFPACK - xGMG]]
Max[Abs[xUMFPACK - xTAUCS]]
Max[Abs[xUMFPACK - xPardiso]]
Max[Abs[xUMFPACK - xKrylov]]

3.90012*10^-10

1.14953*10^-9

2.45955*10^-10

6.27234*10^-10

Ils ont tous une précision comparable. Ainsi, ce simple solveur multigrille, implémenté en un long après-midi, semble être au moins à égalité avec Pardiso. Pas mal, n'est-ce pas ?

Les résolutions multiples sont toujours plus rapides avec les solveurs directs sur les grilles 2D.

Une fois que les factorisations des solveurs directs sont calculées et stockées dans le fichier LinearSolveFunction objets, les résolutions réelles (substitutions avant et arrière du triangle) sont beaucoup plus plus rapides. Cependant, ce n'est pas nécessairement le spectre d'utilisation des méthodes itératives.
Quoi qu'il en soit, voici quelques timings :

solUMFPACK = Quiet[LinearSolve[A]]; // AbsoluteTiming // First
xUMFPACK = solUMFPACK[b]; // AbsoluteTiming // First
solTAUCS = LinearSolve[A, "Method" -> "Cholesky"]; // AbsoluteTiming // First
xTAUCS = solTAUCS[b]; // AbsoluteTiming // First
solPardiso = LinearSolve[A, "Method" -> "Pardiso"]; // AbsoluteTiming // First
xPardiso = solPardiso[b]; // AbsoluteTiming // First

6.07364

0.142823

7.28346

0.183195

2.13817

0.236214

Notez que j'ai utilisé Quiet pour UMFPACK car il se plaint d'un mauvais numéro de condition du système et la gestion des erreurs ajouterait environ 20( !) secondes aux timings. Il n'y a cependant aucun problème avec les erreurs numériques :

Max[Abs[xGMG - xUMFPACK]]
Max[Abs[xGMG - xTAUCS]]
Max[Abs[xGMG - xPardiso]]

3.90012*10^-10

7.59533*10^-10

1.44077*10^-10

Remarques

Le succès des solveurs multigrilles dépend fortement des lisseurs. Ils doivent être bon marché mais aussi efficaces pour se débarrasser des oscillations dans les résidus. Utilisation de CGLinearSolve comme lisseur ne fonctionnera probablement que pour les matrices de système (semi)définies positives. Je pourrais ajouter d'autres lisseurs plus tard.

C'est une implémentation plutôt prématurée et pas entièrement testée. Par exemple, j'aimerais aussi la tester sur des mailles tétraédriques où les méthodes multigrilles sont censées briller. Mais actuellement, je n'ai pas de belles routines pour créer des opérateurs de prolongation.

Exemple 3D

Le problème avec les solveurs directs est qu'à partir de 3 dimensions, leurs performances pour traiter des matrices issues d'EDP chutent rapidement. C'est pourquoi je voulais montrer au moins un exemple en 3 dimensions. Comme il n'y a pas d'analogon immédiat pour la subdivision en boucles des mailles tétraédriques, j'utilise plutôt des mailles hexaédriques.

Préparation

Voici quelques fonctions d'aide pour générer la hiérarchie de la grille et les mappages de prolongation. La prolongation est effectuée par interpolation trilinéaire hexagonale. Cela convient bien puisque nous allons utiliser des mailles du premier ordre.

Needs["NDSolve`FEM`"]
cubemesh[n_] := Module[{R},
  R = ArrayMesh[ConstantArray[1, ConstantArray[n, 3]], DataRange -> ConstantArray[{0, 1}, 3]];
  ToElementMesh[
   "Coordinates" -> MeshCoordinates[R],
   "MeshElements" -> {HexahedronElement[ MeshCells[R, 3, "Multicells" -> True][[1, 1]]]},
   "MeshOrder" -> 1, "NodeReordering" -> False, MeshQualityGoal -> 0
   ]
  ]

getEdges = Compile[{{i, _Integer}, {idx, _Integer, 1}},
   Table[{i, Compile`GetElement[idx, j]}, {j, 1, Length[idx]}],
   CompilationTarget -> "C",
   RuntimeAttributes -> {Listable},
   Parallelization -> True,
   RuntimeOptions -> "Speed"
   ];

getProlongationOperator[Rfine_ElementMesh, Rcoarse_ElementMesh, h_] :=
   Module[{pfine, pcoarse},
   pfine = Rfine["Coordinates"];
   pcoarse = Rcoarse["Coordinates"];
   # SparseArray[1./Total[#, {2}]] &@SparseArray[
     Join @@ getEdges[
        Range[Length[pfine]],
        Nearest[pcoarse -> Automatic, pfine, {∞, h 1.1}, 
         DistanceFunction -> ChessboardDistance]
        ] -> 1.,
     {Length[pfine], Length[pcoarse]}, 0.
     ]
   ];

Ceci crée la hiérarchie de grille réelle et les mappages de prolongation.

dList = Range[6, 2, -1];
nList = 2^dList;
hList = 1./(2^(dList));
RList = cubemesh /@ nList; // AbsoluteTiming // First
J = Table[ getProlongationOperator[RList[[i]], RList[[i + 1]], 
      hList[[i]]], {i, 1, Length[RList] - 1}]; // AbsoluteTiming // First

3.84804

0.603694

Encore une fois, nous résolvons un problème elliptique linéaire avec des conditions aux limites homogènes de Neumann (c'est plus facile à mettre en œuvre que les conditions de Dirichlet). De plus, j'ai pensé que ce serait une bonne idée de prescrire une solution analytique, afin que nous puissions discrétiser son côté droit de l'équation, résoudre l'EDP discrétisée, et comparer avec la solution analytique à la fin. (Notez qu'il est essentiel que v ci-dessous satisfasse les conditions limites homogènes de Neumann).

ϵ = 1.;
Quiet[
  XX = {X[[1]], X[[2]], X[[3]]};
  v = X [Function] Cos[5 Pi X[[1]]] Cos[Pi X[[2]]] Cos[3 Pi X[[3]]];
  Δv = X [Function] Evaluate[Tr[D[v[XX], {XX, 2}]]];
  f = X [Function] Evaluate[ϵ v[XX] - Δv[XX]]
  ];

vd = NDSolve`VariableData[{"DependentVariables", "Space"} -> {{u}, {x, y, z}}];
sd = NDSolve`SolutionData[{"Space"} -> {RList[[1]]}];
cdata = InitializePDECoefficients[vd, sd,
   "DiffusionCoefficients" -> {{-IdentityMatrix[3]}},
   "MassCoefficients" -> {{1}},
   "LoadCoefficients" -> {{f[{x, y, z}]}}
   ];
bcdata = InitializeBoundaryConditions[vd, sd, {DirichletCondition[u[x, y, z] == 0., True]}];
mdata = InitializePDEMethodData[vd, sd];
dpde = DiscretizePDE[cdata, mdata, sd]; // AbsoluteTiming // First

{b, L, damping, M} = dpde["All"];
b = Flatten[Normal[b]];
A = L + ϵ M;

2.21493

La maille la plus fine

consiste en 262144 hexaèdres. La matrice du système a une taille {274625, 274625} et contient 7189057 valeurs non nulles.

Timings

Passons maintenant aux timings. Cette fois, nous voyons que le solveur de gradient conjugué (avec le préconditionneur "ILU0") effectue... beaucoup mieux que les solveurs directs :

xUMFPACK = LinearSolve[A, b]; // AbsoluteTiming // First
xPardiso = LinearSolve[A, b, Method -> "Pardiso"]; // AbsoluteTiming // First
solCG = CGLinearSolve[A, b,
     "Tolerance" -> 1. 10^-6,
     "Preconditioner" -> "ILU0"]; // AbsoluteTiming // First
xCG = solCG["Solution"];

141.175

32.0759

1.70319

Je tiens à préciser qu'une grande partie du temps nécessaire à UMFPACK est due à la gestion de la mémoire de l'OS (je n'ai que 16 Go de RAM installée).

Avec des paramètres légèrement affinés, le solveur multigrille géométrique est encore plus performant :

S = GMGLinearSolve[A, J]; // AbsoluteTiming // First
solGMG = S[b, "Tolerance" -> 1. 10^-4,
     "StartingVectorSmoothingCounts" -> 6,
     "PreSmoothingCounts" -> 4,
     "PostSmoothingCounts" -> 4
     ]; // AbsoluteTiming // First
xGMG = solGMG["Solution"];

0.405304

0.242293

Vous pouvez inspecter [email protected] et voir par vous-même que seulement deux ( !) cycles en V ont été nécessaires. Donc, S peut résoudre la même équation avec des centaines ( !) de côtés droits avant que Pardiso n'ait terminé la factorisation et avant qu'il ne puisse commencer la première résolution réelle. Donc, en pratique, même le LinearSolveFunction générés par LinearSolve compenseront à peine cette différence.

Erreurs

Max[Abs[xUMFPACK - xPardiso]]
Max[Abs[xUMFPACK - xCG]]
Max[Abs[xUMFPACK - xGMG]]

1.42109*10^-14

6.63873*10^-6

3.34638*10^-6

À première vue, cela ne semble pas si bon pour les solveurs itératifs, mais en prenant également en compte l'erreur de discrétisation, ces erreurs sont négligeables :

xTrue = Map[v, RList[[1]]["Coordinates"]];
Max[Abs[xTrue - xPardiso]]/Max[Abs[xTrue]]
Max[Abs[xTrue - xCG]]/Max[Abs[xTrue]]
Max[Abs[xTrue - xGMG]]/Max[Abs[xTrue]]

0.00298862

0.00299525

0.00298605

Vous avez la possibilité de donner de la visibilité à ce post s'il vous a été utile.



Utilisez notre moteur de recherche

Ricerca
Generic filters

Laisser un commentaire

Votre adresse e-mail ne sera pas publiée.