Skip to content

Comment diminuer le temps de calcul lorsque la non-linéarité apparaît dans une équation de conduction thermique en 2D ?

Après avoir observé dans différents référentiels et sites Internet à la fin, nous découvrons la réponse que nous vous apprendrons bientôt.

Solution :

Comme promis, voici mes 6 pence.

Paramètres de base

Needs["NDSolve`FEM`"];
Needs["DifferentialEquations`NDSolveProblems`"];
Needs["DifferentialEquations`NDSolveUtilities`"];

Lr = 2*10^-2;(*dimension of computational domain in r-direction*)
Lz = 10^-2;(*dimension of computational domain in z-direction*)
mesh = ToElementMesh[FullRegion[2], {{0, Lr}, {0, Lz}}, MaxCellMeasure -> {"Length" -> Lr/50}, "MeshOrder" -> 1]
mesh["Wireframe"]

lambda = 22.;         (*heat conductivity*)
density = 7200.;      (*density*)
Cs = 700.;            (*specific heat capacity of solid*) 
Cl = 780.;            (*specific heat capacity of liquid*)      
LatHeat = 272.*10^3;  (*latent heat of fusion*) 
Tliq = 1812.;         (*melting temperature*)
MeltRange = 100.;     (*melting range*)
To = 300.;            (*initial temperature*)       
SPow = 1000.;         (*source power*) 
R = Lr/4.;            (*radius of heat source spot*)
a = Log[100.]/R^2;            
qo = (SPow*a)/Pi; 
q[r_] := qo*Exp[-r^2*a]; (*heat flux distribution*)        
tau = 10^-3;         (*time step size*)
ProcDur = 0.2;       (*process duration*)

Heviside[x_, delta_] := Piecewise[{{0, 
       Abs[x] < -delta}, {0.5*(1 + x/delta + 1/Pi*Sin[(Pi*x)/delta]), 
       Abs[x] <= delta}, {1, x > delta}}];

HevisideDeriv[x_, delta_] := Piecewise[{{0, 
       Abs[x] > delta}, {1/(2*delta)*(1 + Cos[(Pi*x)/delta]), 
       Abs[x] <= delta}}];

EffectHeatCapac[tempr_] := Module[{phase}, 
   phase = Heviside[tempr - Tliq, MeltRange/2];
   Cs*(1 - phase) + Cl*phase + LatHeat*HevisideDeriv[tempr - Tliq, 0.5*MeltRange]];

Versions compilées des fonctions de Heaviside lissées.

cHeaviside = Compile[{{x, _Real}, {delta, _Real}},
   Piecewise[{
     {0., 
      Abs[x] < -delta}, {0.5*(1 + x/delta + 1./Pi*Sin[(Pi*x)/delta]), 
      Abs[x] <= delta}, {1., x > delta}}
    ],
   CompilationTarget -> "C",
   RuntimeAttributes -> {Listable},
   Parallelization -> True
   ];
cHeavisideDeriv = Compile[{{x, _Real}, {delta, _Real}},
   Piecewise[{
     {0., Abs[x] > delta},
     {1./(2*delta)*(1. + Cos[(Pi*x)/delta]), Abs[x] <= delta}}
    ],
   CompilationTarget -> "C",
   RuntimeAttributes -> {Listable},
   Parallelization -> True
   ];
cEffectHeatCapac[tempr_] := 
  With[{phase = cHeaviside[tempr - Tliq, MeltRange/2]},
   Cs*(1 - phase) + Cl*phase + LatHeat*cHeavisideDeriv[tempr - Tliq, 0.5*MeltRange]
   ];

Une routine d'assembleur de matrice rapide

Copié d'ici.

SetAttributes[AssemblyFunction, HoldAll];

Assembly::expected = "Values list has `2` elements. Expected are `1` elements. Returning  prototype.";

Assemble[pat_?MatrixQ, dims_, background_: 0.] := 
  Module[{pa, c, ci, rp, pos}, 
   pa = SparseArray`[email protected][pat -> _, dims];
   rp = pa["RowPointers"];
   ci = pa["ColumnIndices"];
   c = Length[ci];
   pos = cLookupAssemblyPositions[Range[c], rp, Flatten[ci], pat];
   Module[{a},
    a = <|
      "Dimensions" -> dims,
      "Positions" -> pos,
      "RowPointers" -> rp,
      "ColumnIndices" -> ci,
      "Background" -> background,
      "Length" -> c
      |>;
    AssemblyFunction @@ {a}]
   ];

AssemblyFunction /: a_AssemblyFunction[vals0_] := 
  Module[{len, expected, dims, u, vals, dat},
   dat = a[[1]];
   If[VectorQ[vals0], vals = vals0, vals = Flatten[vals0]];
   len = Length[vals];
   expected = Length[dat[["Positions"]]];
   dims = dat[["Dimensions"]];
   If[len === expected, 
    If[Length[dims] == 1, 
     u = ConstantArray[0., dims[[1]]];
     u[[dat[["ColumnIndices"]]]] = AssembleDenseVector[dat[["Positions"]], vals, {dat[["Length"]]}];
     u, 
     SparseArray @@ {Automatic, dims, 
       dat[["Background"]], {1, {dat[["RowPointers"]], 
         dat[["ColumnIndices"]]}, 
        AssembleDenseVector[dat[["Positions"]], 
         vals, {dat[["Length"]]}]}}
     ],
    Message[Assembly::expected, expected, len];
    Abort[]]
   ];

cLookupAssemblyPositions = 
  Compile[{{vals, _Integer, 1}, {rp, _Integer, 1}, {ci, _Integer, 1}, {pat, _Integer, 1}},
   Block[{k, c, i, j},
    i = Compile`GetElement[pat, 1];
    j = Compile`GetElement[pat, 2];
    k = Compile`GetElement[rp, i] + 1;
    c = Compile`GetElement[rp, i + 1];
    While[k < c + 1 && Compile`GetElement[ci, k] != j,
     ++k
     ];
    Compile`GetElement[vals, k]
    ],
   RuntimeAttributes -> {Listable},
   Parallelization -> True,
   CompilationTarget -> "C",
   RuntimeOptions -> "Speed"
   ];

AssembleDenseVector = 
  Compile[{{ilist, _Integer, 1}, {values, _Real, 1}, {dims, _Integer, 1}}, Block[{A}, A = Table[0., {Compile`GetElement[dims, 1]}];
    Do[A[[Compile`GetElement[ilist, i]]] += 
      Compile`GetElement[values, i], {i, 1, Length[values]}];
    A
    ],
   CompilationTarget -> "C",
   RuntimeOptions -> "Speed"
   ];

Code d'assemblage de matrice d'amortissement

Principalement de la rétro-ingénierie, donc je ne suis effectivement pas sûr à 100% que cela fait ce qu'il faut....

Pour autant que je l'ai obtenu, la matrice d'amortissement par rapport à la fonction... $f colon varOmega to mathbb{R}$ devrait coder la forme bilinéaire

$$(u,v) mapsto int_{varOmega} u(x) , v(x) , f(x) , mathrm{d} x.$$ en termes de fonctions de base FEM. Comme les fonctions de base FEM ont un support très local, on passe par les éléments finis du maillage (des quads dans notre cas) et on calcule les contributions locales à la forme bilinéaire globale. Ensuite, il s'agit de jongler avec les indices pour assembler la forme bilinéaire.

Cela suppose une interpolation bi-linéaire sur les quads et emploie la quadrature gaussienne avec 2 points d'intégration par dimension pour l'intégration. (Pour les maillages triangulaires ou tétraédriques, l'intégration exacte peut être utilisée à la place).

(* for each quad, `getWeakLaplaceCombinatoricsQuad` is supposed to produce the $i-j$-indices of each of the 16 entries of the local $4 times 4$ metrix within the global matrix.*)
getWeakLaplaceCombinatoricsQuad = Block[{q},
   With[{code = Flatten[Table[Table[{
          Compile`GetElement[q, i],
          Compile`GetElement[q, j]
          }, {i, 1, 4}], {j, 1, 4}], 1]},
    Compile[{{q, _Integer, 1}},
     code,
     CompilationTarget -> "C",
     RuntimeAttributes -> {Listable},
     Parallelization -> True,
     RuntimeOptions -> "Speed"
     ]
    ]
   ];

(* this snippet computes the symbolic expression for the local matrices and then compiles it into the function `getLocalDampingMatrices`*) 
Block[{dim, PP, UU, FF, p, u, f, integrant, x, ω, localmatrix},
  dim = 2;
  PP = Table[Compile`GetElement[P, i, j], {i, 1, 4}, {j, 1, dim}];
  UU = Table[Compile`GetElement[U, i], {i, 1, 4}];
  FF = Table[Compile`GetElement[F, i], {i, 1, 4}];

  (* bi-linear interpolation of the quadrilateral; maps the standard quare onto the quadrilateral defined by PP[[1]], PP[[2]], PP[[3]], PP[[3]]*)
  p = {s, t} [Function] (PP[[1]] (1 - s) + s PP[[2]]) (1 - t) + t (PP[[4]] (1 - s) + s PP[[3]]);
  (* bi-linear interpolation of the function values of u*)
  u = {s, t} [Function] (UU[[1]] (1 - s) + s UU[[2]]) (1 - t) + t (UU[[4]] (1 - s) + s UU[[3]]);
  (* bi-linear interpolation of the function values of f*)
  f = {s, t} [Function] (FF[[1]] (1 - s) + s FF[[2]]) (1 - t) + t integrant = {s, t} [Function] Evaluate[f[s, t] u[s, t]^2 Abs[Det[D[p[s, t], {{s, t}, 1}]]]];
  {x, ω} = Most[NIntegrate`GaussRuleData[2, MachinePrecision]];

  (* using `D` to obtain the local matrix from its quadratic form*)
  localmatrix = 1/2 D[
     Flatten[KroneckerProduct[ω, ω]].integrant @@@ Tuples[x, 2],
     {UU, 2}
     ];

  (* `getLocalDampingMatrices` computes the local $4 times 4$-matrices from the quad vertex coordinates `P` (supposed to be a $4 times 2$-matrix) and from the function values `F` (supposed to be a $4$-vector) *) 
  getLocalDampingMatrices = With[{code = localmatrix},
    Compile[{{P, _Real, 2}, {F, _Real, 1}},
     code,
     CompilationTarget -> "C",
     RuntimeAttributes -> {Listable},
     Parallelization -> True,
     RuntimeOptions -> "Speed"
     ]
    ];
  ];

getDampingMatrix[assembler_AssemblyFunction, quads_, quaddata_, fvals_] := 
  Module[{fdata, localmatrices},
   fdata = Partition[fvals[[Flatten[quads]]], 4];
   localmatrices = getLocalDampingMatrices[quaddata, fdata];
   assembler[Flatten[localmatrices]]
   ];

La fonction getDampingMatrix mange un AssemblyFunction objet assembler_, la liste quads de tous les quads (sous la forme d'une liste de 4-vecteurs des indices des sommets), la liste quaddata (une liste de 4$ times 2-matrice avec les positions des sommets, et une liste de fvals avec les valeurs de la fonction $f$ aux sommets de la maille. Il crache la matrice d'amortissement complètement assemblée.

Utilisation de DiscretizePDE une seule fois

Cela nécessite l'ancienne implémentation de EffectHeatCapac.

u =.
vd = NDSolve`VariableData[{"DependentVariables" -> {u}, "Space" -> {r, z}, "Time" -> t}];
sd = NDSolve`SolutionData[{"Space", "Time"} -> {ToNumericalRegion[mesh], 0.}];

DirichCond = DirichletCondition[u[t, r, z] == To, z == 0];
NeumCond = NeumannValue[q[r], z == Lz];
initBCs = InitializeBoundaryConditions[vd, sd, {{DirichCond, NeumCond}}];
methodData = InitializePDEMethodData[vd, sd];
discreteBCs = DiscretizeBoundaryConditions[initBCs, methodData, sd];

x0 = ConstantArray[To, {methodData["DegreesOfFreedom"]}];
TemprField = ElementMeshInterpolation[{mesh}, x0];
NumTimeStep = Floor[ProcDur/tau];

pdeCoefficients = InitializePDECoefficients[vd, sd,
   "ConvectionCoefficients" -> {{{{-(lambda/r), 0}}}},
   "DiffusionCoefficients" -> {{-lambda*IdentityMatrix[2]}}, 
   "DampingCoefficients" -> {{EffectHeatCapac[TemprField[r, z]] density}}
   ];
discretePDE = DiscretizePDE[pdeCoefficients, methodData, sd];
{load, stiffness, damping, mass} = discretePDE["SystemMatrices"];
DeployBoundaryConditions[{load, stiffness, damping}, discreteBCs];

Exécution de la simulation

En supprimant les goulots d'étranglement DiscretizePDE et (beaucoup plus sévèrement) ElementMeshInterpolationla boucle ne nécessite plus que 0,32 seconde pour s'exécuter. Nous profitons également du fait que, en utilisant la fonction AssemblyFunctionassemblernous n'avons pas besoin de recalculer les motifs de tableaux épars. De plus, l'utilisation d'une syntaxe non documentée pour la fonction SparseArray permet de contourner certaines dégradations supplémentaires des performances.

Donc, c'est maintenant plus rapide par un facteur de 100.

x = x0;
taustiffness = tau stiffness;
tauload = tau Flatten[load];

quads = mesh["MeshElements"][[1, 1]];
quaddata = Partition[mesh["Coordinates"][[Flatten[quads]]], 4];
assembler = Assemble[Flatten[getWeakLaplaceCombinatoricsQuad[quads], 1], {1, 1} Length[mesh["Coordinates"]]];

Do[
    damping = getDampingMatrix[assembler, quads, quaddata, cEffectHeatCapac[x] density];
    DeployBoundaryConditions[{load, stiffness, damping}, discreteBCs];
    A = damping + taustiffness;
    b = tauload + damping.x;
    x = LinearSolve[A, b, Method -> {"Krylov",
        Method -> "BiCGSTAB",
        "Preconditioner" -> "ILU0",
        "StartingVector" -> x
        }
      ];
    ,
    {i, 1, NumTimeStep}]; // AbsoluteTiming // First

0.325719

Utilisation de ElementMeshInterpolation seulement une fois à la fin pour le tracé

TemprField = ElementMeshInterpolation[{mesh}, x];

ContourPlot[TemprField[r, z], {r, z} ∈ mesh,
 AspectRatio -> Lz/Lr,
 ColorFunction -> "TemperatureMap",
 Contours -> 50,
 PlotRange -> All,
 PlotLegends -> Placed[Automatic, After],
 FrameLabel -> {"r", "z"},
 PlotPoints -> 50,
 PlotLabel -> "Temperature field",
 BaseStyle -> 16]

Entrez la description de l'image ici

Addendum

Après l'exécution

fvals = cEffectHeatCapac[x] density;
fdata = Partition[fvals[[Flatten[quads]]], 4];
localmatrices = getLocalDampingMatrices[quaddata, fdata];

la ligne

assembler[localmatrices];

est fondamentalement équivalente à l'utilisation de SparseArray pour l'assemblage additif, comme suit :

(* switching to additive matrix assembly *)
SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> Total}];
pat = Join @@ getWeakLaplaceCombinatoricsQuad[quads];
SparseArray[pat -> Flatten[localmatrices], {1, 1} Length[fvals], 0.];

Peut-être que cela aide à comprendre comment getWeakLaplaceCombinatoricsQuad et getLocalDampingMatrices sont liés.

Addendum II

J'ai implémenté une interface un peu plus lisse pour les maillages simpliciaux de dimensions arbitraires ici.

Supposons donc que nous ayons commencé avec le maillage triangulaire suivant :

mesh = ToElementMesh[FullRegion[2], {{0, Lr}, {0, Lz}}, 
   MaxCellMeasure -> {"Length" -> Lr/50}, "MeshOrder" -> 1,
   MeshElementType -> TriangleElement];

Ensuite, on doit convertir le maillage une fois en un MeshRegion.

Ω = MeshRegion[mesh];

et au lieu de

damping = getDampingMatrix[assembler, quads, quaddata, cEffectHeatCapac[x] density];

ainsi que la définition de assembler, quads, quaddataetc., on peut simplement utiliser

damping = RegionReactionMatrix[Ω, cEffectHeatCapac[x] density]

dans les Do-boucle.

J'ai réussi à réduire le temps par 2,5 fois + j'ai ajouté la possibilité d'afficher la température en fonction de l'heure. J'ai utilisé Do[] au lieu de For[] et Interpolation[] au lieu de Module[]. Nous pouvons encore accélérer le code.

Needs["NDSolve`FEM`"];
Needs["DifferentialEquations`NDSolveProblems`"];
Needs["DifferentialEquations`NDSolveUtilities`"];
Lr = 2*10^-2;(*dimension of computational domain in r-direction*)Lz = 
 10^-2;(*dimension of computational domain in z-direction*)mesh = 
 ToElementMesh[FullRegion[2], {{0, Lr}, {0, Lz}}, 
  MaxCellMeasure -> {"Length" -> Lr/50}, "MeshOrder" -> 1]
mesh["Wireframe"]
lambda = 22;(*heat conductivity*)density = 7200;(*density*)Cs = 
700;(*specific heat capacity of solid*)Cl = 780;(*specific heat 
capacity of liquid*)LatHeat = 
 272*10^3;(*latent heat of fusion*)Tliq = 1812;(*melting 
temperature*)MeltRange = 100;(*melting range*)To = 300;(*initial 
temperature*)SPow = 1000;(*source power*)R = 
 Lr/4;(*radius of heat source spot*)a = Log[100]/R^2;
qo = (SPow*a)/Pi;
q[r_] := qo*Exp[-r^2*a];(*heat flux distribution*)tau = 
 10^-3;(*time step size*)ProcDur = 0.2;(*process duration*)
Heviside[x_, delta_] := 
 Module[{res}, 
  res = Piecewise[{{0, 
      Abs[x] < -delta}, {0.5*(1 + x/delta + 1/Pi*Sin[(Pi*x)/delta]), 
      Abs[x] <= delta}, {1, x > delta}}];
  res]
HevisideDeriv[x_, delta_] := 
 Module[{res}, 
  res = Piecewise[{{0, 
      Abs[x] > delta}, {1/(2*delta)*(1 + Cos[(Pi*x)/delta]), 
      Abs[x] <= delta}}];
  res]
EffectHeatCapac[tempr_] := 
 Module[{phase}, phase = Heviside[tempr - Tliq, MeltRange/2];
  Cs*(1 - phase) + Cl*phase + 
   LatHeat*HevisideDeriv[tempr - Tliq, 0.5*MeltRange]]
ehc = Interpolation[
   Table[{x, EffectHeatCapac[x]}, {x, To - 100, 4000, 1}]];
ts = AbsoluteTime[];

NumTimeStep = Floor[ProcDur/tau];

vd = NDSolve`VariableData[{"DependentVariables" -> {u}, 
    "Space" -> {r, z}, "Time" -> t}];
sd = NDSolve`SolutionData[{"Space", 
     "Time"} -> {ToNumericalRegion[mesh], 0.}];

DirichCond = DirichletCondition[u[t, r, z] == To, z == 0];
NeumCond = NeumannValue[q[r], z == Lz];
initBCs = 
  InitializeBoundaryConditions[vd, sd, {{DirichCond, NeumCond}}];
methodData = InitializePDEMethodData[vd, sd];
discreteBCs = DiscretizeBoundaryConditions[initBCs, methodData, sd];
xlast = Table[{To}, {methodData["DegreesOfFreedom"]}];
TemprField[0] = ElementMeshInterpolation[{mesh}, xlast];
Do[(*(*Setting of PDE coefficients for linear 
problem*)pdeCoefficients=InitializePDECoefficients[vd,sd,
"ConvectionCoefficients"[Rule]{{{{-lambda/r,0}}}},
"DiffusionCoefficients"[Rule]{{-lambda*IdentityMatrix[2]}},
"DampingCoefficients"[Rule]{{Cs*density}}];*)(*Setting of PDE 
coefficients for nonlinear problem*)
 pdeCoefficients = 
  InitializePDECoefficients[vd, sd, 
   "ConvectionCoefficients" -> {{{{-(lambda/r), 0}}}}, 
   "DiffusionCoefficients" -> {{-lambda*IdentityMatrix[2]}}, 
   "DampingCoefficients" -> {{ehc[TemprField[i - 1][r, z]]*density}}];
 discretePDE = DiscretizePDE[pdeCoefficients, methodData, sd];
 {load, stiffness, damping, mass} = discretePDE["SystemMatrices"];
 DeployBoundaryConditions[{load, stiffness, damping}, discreteBCs];
 A = damping/tau + stiffness;
 b = load + damping.xlast/tau;
 x = LinearSolve[A, b, 
   Method -> {"Krylov", Method -> "BiCGSTAB", 
     "Preconditioner" -> "ILU0", 
     "StartingVector" -> Flatten[xlast, 1]}];
 TemprField[i] = ElementMeshInterpolation[{mesh}, x];
 xlast = x;, {i, 1, NumTimeStep}]
te = AbsoluteTime[];
te - ts

Voici le temps pour l'ancien et le nouveau code 39.4973561 et 15.4960282 respectivement (sur mon ASUS ZenBook).Pour réduire encore le temps, utilisez l'option MeshRefinementFunction:

f = Function[{vertices, area}, 
  Block[{r, z}, {r, z} = Mean[vertices]; 
   If[r^2 + (z - Lz)^2 <= (Lr/4)^2, area > (Lr/50)^2, 
    area > (Lr/
        15)^2]]];
mesh = 
 ToElementMesh[FullRegion[2], {{0, Lr}, {0, Lz}}, "MeshOrder" -> 1, 
  MeshRefinementFunction -> f]
mesh["Wireframe"]

Pour cette option le temps est de 8.8878213

{ContourPlot[TemprField[NumTimeStep][r, z], {r, 0, Lr}, {z, 0, Lz}, 
  PlotRange -> All, ColorFunction -> "TemperatureMap", 
  PlotLegends -> Automatic, FrameLabel -> Automatic], 
 ListPlot[Table[{tau*i, TemprField[i][.001, Lz]}, {i, 0, 
    NumTimeStep}], AxesLabel -> {"t", "T"}]}

Fig. 1
Grâce à Henrik Schumacher, nous pouvons encore accélérer le code. J'ai légèrement modifié son code dans le cas de l'utilisation du "WVM" et pour afficher le champ de température à chaque étape.

Needs["NDSolve`FEM`"];
Needs["DifferentialEquations`NDSolveProblems`"];
Needs["DifferentialEquations`NDSolveUtilities`"];

SetAttributes[AssemblyFunction, HoldAll];

Assembly::expected = 
  "Values list has `2` elements. Expected are `1` elements. Returning 
 prototype.";

Assemble[pat_?MatrixQ, dims_, background_: 0.] := 
  Module[{pa, c, ci, rp, pos}, 
   pa = SparseArray`SparseArraySor[email protected][pat -> _, dims];
   rp = pa["RowPointers"];
   ci = pa["ColumnIndices"];
   c = Length[ci];
   pos = cLookupAssemblyPositions[Range[c], rp, Flatten[ci], pat];
   Module[{a}, 
    a = <|"Dimensions" -> dims, "Positions" -> pos, 
      "RowPointers" -> rp, "ColumnIndices" -> ci, 
      "Background" -> background, "Length" -> c|>;
    AssemblyFunction @@ {a}]];

AssemblyFunction /: a_AssemblyFunction[vals0_] := 
  Module[{len, expected, dims, u, vals, dat}, dat = a[[1]];
   If[VectorQ[vals0], vals = vals0, vals = Flatten[vals0]];
   len = Length[vals];
   expected = Length[dat[["Positions"]]];
   dims = dat[["Dimensions"]];
   If[len === expected, 
    If[Length[dims] == 1, u = ConstantArray[0., dims[[1]]];
     u[[dat[["ColumnIndices"]]]] = 
      AssembleDenseVector[dat[["Positions"]], vals, {dat[["Length"]]}];
     u, SparseArray @@ {Automatic, dims, 
       dat[["Background"]], {1, {dat[["RowPointers"]], 
         dat[["ColumnIndices"]]}, 
        AssembleDenseVector[dat[["Positions"]], 
         vals, {dat[["Length"]]}]}}], 
    Message[Assembly::expected, expected, len];
    Abort[]]];

cLookupAssemblyPositions = 
  Compile[{{vals, _Integer, 1}, {rp, _Integer, 1}, {ci, _Integer, 
     1}, {pat, _Integer, 1}}, 
   Block[{k, c, i, j}, i = Compile`GetElement[pat, 1];
    j = Compile`GetElement[pat, 2];
    k = Compile`GetElement[rp, i] + 1;
    c = Compile`GetElement[rp, i + 1];
    While[k < c + 1 && Compile`GetElement[ci, k] != j, ++k];
    Compile`GetElement[vals, k]], RuntimeAttributes -> {Listable}, 
   Parallelization -> True, CompilationTarget -> "WVM", 
   RuntimeOptions -> "Speed"];

AssembleDenseVector = 
  Compile[{{ilist, _Integer, 1}, {values, _Real, 1}, {dims, _Integer, 
     1}}, Block[{A}, A = Table[0., {Compile`GetElement[dims, 1]}];
    Do[A[[Compile`GetElement[ilist, i]]] += 
      Compile`GetElement[values, i], {i, 1, Length[values]}];
    A], CompilationTarget -> "WVM", RuntimeOptions -> "Speed"];
getWeakLaplaceCombinatoricsQuad =   
Block[{q}, 
   With[{code = 
      Flatten[Table[
        Table[{Compile`GetElement[q, i], 
          Compile`GetElement[q, j]}, {i, 1, 4}], {j, 1, 4}], 1]}, 
    Compile[{{q, _Integer, 1}}, code, CompilationTarget -> "WVM", 
     RuntimeAttributes -> {Listable}, Parallelization -> True, 
     RuntimeOptions -> "Speed"]]];

Block[{dim, PP, UU, FF, p, u, f, integrant, x, [Omega], localmatrix},
   dim = 2;
  PP = Table[Compile`GetElement[P, i, j], {i, 1, 4}, {j, 1, dim}];
  UU = Table[Compile`GetElement[U, i], {i, 1, 4}];
  FF = Table[Compile`GetElement[F, i], {i, 1, 4}];
  p = {s, t} [Function] (PP[[1]] (1 - s) + s PP[[2]]) (1 - t) + 
     t (PP[[4]] (1 - s) + s PP[[3]]);
  u = {s, t} [Function] (UU[[1]] (1 - s) + s UU[[2]]) (1 - t) + 
     t (UU[[4]] (1 - s) + s UU[[3]]);
  f = {s, t} [Function] (FF[[1]] (1 - s) + s FF[[2]]) (1 - t) + 
     t (FF[[4]] (1 - s) + s FF[[3]]);
  integrant = {s, t} [Function] 
    Evaluate[f[s, t] u[s, t]^2 Abs[Det[D[p[s, t], {{s, t}, 1}]]]];
  {x, [Omega]} = Most[NIntegrate`GaussRuleData[2, MachinePrecision]];
  localmatrix = 
   1/2 D[Flatten[KroneckerProduct[[Omega], [Omega]]].integrant @@@ 
       Tuples[x, 2], {UU, 2}];
  getLocalDampingMatrices = 
   With[{code = localmatrix}, 
    Compile[{{P, _Real, 2}, {F, _Real, 1}}, code, 
     CompilationTarget -> "WVM", RuntimeAttributes -> {Listable}, 
     Parallelization -> True, RuntimeOptions -> "Speed"]];];

getDampingMatrix[assembler_, quads_, quaddata_, vals_] := 
  Module[{fvals, fdata, localmatrices}, 
   fvals = cEffectHeatCapac[Flatten[vals]]*density;
   fdata = Partition[fvals[[Flatten[quads]]], 4];
   localmatrices = getLocalDampingMatrices[quaddata, fdata];
   assembler[Flatten[localmatrices]]];
Lr = 2*10^-2;(*dimension of computational domain in r-direction*)Lz = 
 10^-2;(*dimension of computational domain in z-direction*)mesh = 
 ToElementMesh[FullRegion[2], {{0, Lr}, {0, Lz}}, 
  MaxCellMeasure -> {"Length" -> Lr/50}, "MeshOrder" -> 1]
mesh["Wireframe"]

lambda = 22.;(*heat conductivity*)density = 7200.;(*density*)Cs = 
700.;(*specific heat capacity of solid*)Cl = 780.;(*specific heat 
capacity of liquid*)LatHeat = 
 272.*10^3;(*latent heat of fusion*)Tliq = 1812.;(*melting 
temperature*)MeltRange = 100.;(*melting range*)To = 300.;(*initial 
temperature*)SPow = 1000.;(*source power*)R = 
 Lr/4.;(*radius of heat source spot*)a = Log[100.]/R^2;
qo = (SPow*a)/Pi;
q[r_] := qo*Exp[-r^2*a];(*heat flux distribution*)tau = 
 10^-3;(*time step size*)ProcDur = 0.2;(*process duration*)
Heviside[x_, delta_] := 
 Piecewise[{{0, 
    Abs[x] < -delta}, {0.5*(1 + x/delta + 1/Pi*Sin[(Pi*x)/delta]), 
    Abs[x] <= delta}, {1, x > delta}}];

HevisideDeriv[x_, delta_] := 
  Piecewise[{{0, 
     Abs[x] > delta}, {1/(2*delta)*(1 + Cos[(Pi*x)/delta]), 
     Abs[x] <= delta}}];

EffectHeatCapac[tempr_] := 
  Module[{phase}, phase = Heviside[tempr - Tliq, MeltRange/2];
   Cs*(1 - phase) + Cl*phase + 
    LatHeat*HevisideDeriv[tempr - Tliq, 0.5*MeltRange]];
cHeaviside = 
  Compile[{{x, _Real}, {delta, _Real}}, 
   Piecewise[{{0., 
      Abs[x] < -delta}, {0.5*(1 + x/delta + 1./Pi*Sin[(Pi*x)/delta]), 
      Abs[x] <= delta}, {1., x > delta}}], CompilationTarget -> "WVM",
    RuntimeAttributes -> {Listable}, Parallelization -> True];
cHeavisideDeriv = 
  Compile[{{x, _Real}, {delta, _Real}}, 
   Piecewise[{{0., 
      Abs[x] > delta}, {1./(2*delta)*(1. + Cos[(Pi*x)/delta]), 
      Abs[x] <= delta}}], CompilationTarget -> "WVM", 
   RuntimeAttributes -> {Listable}, Parallelization -> True];
cEffectHeatCapac[tempr_] := 
  With[{phase = cHeaviside[tempr - Tliq, MeltRange/2]}, 
   Cs*(1 - phase) + Cl*phase + 
    LatHeat*cHeavisideDeriv[tempr - Tliq, 0.5*MeltRange]];
u =.
vd = NDSolve`VariableData[{"DependentVariables" -> {u}, 
    "Space" -> {r, z}, "Time" -> t}];
sd = NDSolve`SolutionData[{"Space", 
     "Time"} -> {ToNumericalRegion[mesh], 0.}];

DirichCond = DirichletCondition[u[t, r, z] == To, z == 0];
NeumCond = NeumannValue[q[r], z == Lz];
initBCs = 
  InitializeBoundaryConditions[vd, sd, {{DirichCond, NeumCond}}];
methodData = InitializePDEMethodData[vd, sd];
discreteBCs = DiscretizeBoundaryConditions[initBCs, methodData, sd];

x0 = ConstantArray[To, {methodData["DegreesOfFreedom"]}];
TemprField = ElementMeshInterpolation[{mesh}, x0];
NumTimeStep = Floor[ProcDur/tau];

pdeCoefficients = 
  InitializePDECoefficients[vd, sd, 
   "ConvectionCoefficients" -> {{{{-(lambda/r), 0}}}}, 
   "DiffusionCoefficients" -> {{-lambda*IdentityMatrix[2]}}, 
   "DampingCoefficients" -> {{EffectHeatCapac[
        TemprField[r, z]] density}}];
discretePDE = DiscretizePDE[pdeCoefficients, methodData, sd];
{load, stiffness, damping, mass} = discretePDE["SystemMatrices"];
DeployBoundaryConditions[{load, stiffness, damping}, discreteBCs];
x = x0;
X[0] = x;
taustiffness = tau stiffness;
tauload = tau Flatten[load];

quads = mesh["MeshElements"][[1, 1]];
quaddata = Partition[mesh["Coordinates"][[Flatten[quads]]], 4];
assembler = 
  Assemble[Flatten[getWeakLaplaceCombinatoricsQuad[quads], 
    1], {1, 1} Length[mesh["Coordinates"]]];

Do[damping = getDampingMatrix[assembler, quads, quaddata, x];
    DeployBoundaryConditions[{load, stiffness, damping}, discreteBCs];
    A = damping + taustiffness;
    b = tauload + damping.x;
    x = LinearSolve[A, b, 
      Method -> {"Krylov", Method -> "BiCGSTAB", 
        "Preconditioner" -> "ILU0", "StartingVector" -> x, 
        "Tolerance" -> 0.00001}]; X[i] = x;, {i, 1, NumTimeStep}]; // 
  AbsoluteTiming // First

Ici nous avons le temps 0.723424 et la température à chaque étape

T[i_] := ElementMeshInterpolation[{mesh}, X[i]]

ContourPlot[T[NumTimeStep][r, z], {r, z} [Element] mesh, 
 AspectRatio -> Lz/Lr, ColorFunction -> "TemperatureMap", 
 PlotLegends -> Automatic, PlotRange -> All, Contours -> 20]

ListPlot[Table[{i*tau, T[i][.001, Lz]}, {i, 0, NumTimeStep}]]

Fig2

Section des critiques et des évaluations



Utilisez notre moteur de recherche

Ricerca
Generic filters

Laisser un commentaire

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