FreeFOAM The Cross-Platform CFD Toolkit
MULESTemplates.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 1991-2011 OpenCFD Ltd.
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "MULES.H"
27 #include <finiteVolume/upwind.H>
35 #include <OpenFOAM/syncTools.H>
36 
37 #include <finiteVolume/fvCFD.H>
38 
39 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
40 
41 template<class RhoType, class SpType, class SuType>
43 (
44  const RhoType& rho,
46  const surfaceScalarField& phi,
47  surfaceScalarField& phiPsi,
48  const SpType& Sp,
49  const SuType& Su,
50  const scalar psiMax,
51  const scalar psiMin
52 )
53 {
54  Info<< "MULES: Solving for " << psi.name() << endl;
55 
56  const fvMesh& mesh = psi.mesh();
58 
59  surfaceScalarField phiBD = upwind<scalar>(psi.mesh(), phi).flux(psi);
60 
61  surfaceScalarField& phiCorr = phiPsi;
62  phiCorr -= phiBD;
63 
64  scalarField allLambda(mesh.nFaces(), 1.0);
65 
67  (
68  IOobject
69  (
70  "lambda",
71  mesh.time().timeName(),
72  mesh,
73  IOobject::NO_READ,
74  IOobject::NO_WRITE,
75  false
76  ),
77  mesh,
78  dimless,
79  allLambda,
80  false // Use slices for the couples
81  );
82 
83  limiter
84  (
85  allLambda,
86  rho,
87  psi,
88  phiBD,
89  phiCorr,
90  Sp.field(),
91  Su.field(),
92  psiMax,
93  psiMin,
94  3
95  );
96 
97  phiPsi = phiBD + lambda*phiCorr;
98 
99  scalarField& psiIf = psi;
100  const scalarField& psi0 = psi.oldTime();
101  const scalar deltaT = mesh.time().deltaT().value();
102 
103  psiIf = 0.0;
104  fvc::surfaceIntegrate(psiIf, phiPsi);
105 
106  if (mesh.moving())
107  {
108  psiIf =
109  (
110  mesh.Vsc0()*rho.oldTime()*psi0/(deltaT*mesh.Vsc())
111  + Su.field()
112  - psiIf
113  )/(rho/deltaT - Sp.field());
114  }
115  else
116  {
117  psiIf =
118  (
119  rho.oldTime()*psi0/deltaT
120  + Su.field()
121  - psiIf
122  )/(rho/deltaT - Sp.field());
123  }
124 
126 }
127 
128 
129 template<class RhoType, class SpType, class SuType>
131 (
132  const RhoType& rho,
133  volScalarField& psi,
134  const surfaceScalarField& phi,
135  surfaceScalarField& phiPsi,
136  const SpType& Sp,
137  const SuType& Su,
138  const scalar psiMax,
139  const scalar psiMin
140 )
141 {
142  const fvMesh& mesh = psi.mesh();
143 
144  const dictionary& MULEScontrols = mesh.solverDict(psi.name());
145 
146  label maxIter
147  (
148  readLabel(MULEScontrols.lookup("maxIter"))
149  );
150 
151  label nLimiterIter
152  (
153  readLabel(MULEScontrols.lookup("nLimiterIter"))
154  );
155 
156  scalar maxUnboundedness
157  (
158  readScalar(MULEScontrols.lookup("maxUnboundedness"))
159  );
160 
161  scalar CoCoeff
162  (
163  readScalar(MULEScontrols.lookup("CoCoeff"))
164  );
165 
166  scalarField allCoLambda(mesh.nFaces());
167 
168  {
169  surfaceScalarField Cof =
170  mesh.time().deltaT()*mesh.surfaceInterpolation::deltaCoeffs()
171  *mag(phi)/mesh.magSf();
172 
173  slicedSurfaceScalarField CoLambda
174  (
175  IOobject
176  (
177  "CoLambda",
178  mesh.time().timeName(),
179  mesh,
180  IOobject::NO_READ,
181  IOobject::NO_WRITE,
182  false
183  ),
184  mesh,
185  dimless,
186  allCoLambda,
187  false // Use slices for the couples
188  );
189 
190  CoLambda == 1.0/max(CoCoeff*Cof, scalar(1));
191  }
192 
193  scalarField allLambda(allCoLambda);
194  //scalarField allLambda(mesh.nFaces(), 1.0);
195 
197  (
198  IOobject
199  (
200  "lambda",
201  mesh.time().timeName(),
202  mesh,
203  IOobject::NO_READ,
204  IOobject::NO_WRITE,
205  false
206  ),
207  mesh,
208  dimless,
209  allLambda,
210  false // Use slices for the couples
211  );
212 
213  linear<scalar> CDs(mesh);
214  upwind<scalar> UDs(mesh, phi);
215  //fv::uncorrectedSnGrad<scalar> snGrads(mesh);
216 
217  fvScalarMatrix psiConvectionDiffusion
218  (
219  fvm::ddt(rho, psi)
220  + fv::gaussConvectionScheme<scalar>(mesh, phi, UDs).fvmDiv(phi, psi)
221  //- fv::gaussLaplacianScheme<scalar, scalar>(mesh, CDs, snGrads)
222  //.fvmLaplacian(Dpsif, psi)
223  - fvm::Sp(Sp, psi)
224  - Su
225  );
226 
227  surfaceScalarField phiBD = psiConvectionDiffusion.flux();
228 
229  surfaceScalarField& phiCorr = phiPsi;
230  phiCorr -= phiBD;
231 
232  for (label i=0; i<maxIter; i++)
233  {
234  if (i != 0 && i < 4)
235  {
236  allLambda = allCoLambda;
237  }
238 
239  limiter
240  (
241  allLambda,
242  rho,
243  psi,
244  phiBD,
245  phiCorr,
246  Sp.field(),
247  Su.field(),
248  psiMax,
249  psiMin,
250  nLimiterIter
251  );
252 
253  solve
254  (
255  psiConvectionDiffusion + fvc::div(lambda*phiCorr),
256  MULEScontrols
257  );
258 
259  scalar maxPsiM1 = gMax(psi.internalField()) - 1.0;
260  scalar minPsi = gMin(psi.internalField());
261 
262  scalar unboundedness = max(max(maxPsiM1, 0.0), -min(minPsi, 0.0));
263 
264  if (unboundedness < maxUnboundedness)
265  {
266  break;
267  }
268  else
269  {
270  Info<< "MULES: max(" << psi.name() << " - 1) = " << maxPsiM1
271  << " min(" << psi.name() << ") = " << minPsi << endl;
272 
273  phiBD = psiConvectionDiffusion.flux();
274 
275  /*
276  word gammaScheme("div(phi,gamma)");
277  word gammarScheme("div(phirb,gamma)");
278 
279  const surfaceScalarField& phir =
280  mesh.lookupObject<surfaceScalarField>("phir");
281 
282  phiCorr =
283  fvc::flux
284  (
285  phi,
286  psi,
287  gammaScheme
288  )
289  + fvc::flux
290  (
291  -fvc::flux(-phir, scalar(1) - psi, gammarScheme),
292  psi,
293  gammarScheme
294  )
295  - phiBD;
296  */
297  }
298  }
299 
300  phiPsi = psiConvectionDiffusion.flux() + lambda*phiCorr;
301 }
302 
303 
304 template<class RhoType, class SpType, class SuType>
306 (
307  scalarField& allLambda,
308  const RhoType& rho,
309  const volScalarField& psi,
310  const surfaceScalarField& phiBD,
311  const surfaceScalarField& phiCorr,
312  const SpType& Sp,
313  const SuType& Su,
314  const scalar psiMax,
315  const scalar psiMin,
316  const label nLimiterIter
317 )
318 {
319  const scalarField& psiIf = psi;
320  const volScalarField::GeometricBoundaryField& psiBf = psi.boundaryField();
321 
322  const scalarField& psi0 = psi.oldTime();
323 
324  const fvMesh& mesh = psi.mesh();
325 
326  const unallocLabelList& owner = mesh.owner();
327  const unallocLabelList& neighb = mesh.neighbour();
329  const scalarField& V = tVsc();
330  const scalar deltaT = mesh.time().deltaT().value();
331 
332  const scalarField& phiBDIf = phiBD;
333  const surfaceScalarField::GeometricBoundaryField& phiBDBf =
334  phiBD.boundaryField();
335 
336  const scalarField& phiCorrIf = phiCorr;
337  const surfaceScalarField::GeometricBoundaryField& phiCorrBf =
338  phiCorr.boundaryField();
339 
341  (
342  IOobject
343  (
344  "lambda",
345  mesh.time().timeName(),
346  mesh,
347  IOobject::NO_READ,
348  IOobject::NO_WRITE,
349  false
350  ),
351  mesh,
352  dimless,
353  allLambda,
354  false // Use slices for the couples
355  );
356 
357  scalarField& lambdaIf = lambda;
358  surfaceScalarField::GeometricBoundaryField& lambdaBf =
359  lambda.boundaryField();
360 
361  scalarField psiMaxn(psiIf.size(), psiMin);
362  scalarField psiMinn(psiIf.size(), psiMax);
363 
364  scalarField sumPhiBD(psiIf.size(), 0.0);
365 
366  scalarField sumPhip(psiIf.size(), VSMALL);
367  scalarField mSumPhim(psiIf.size(), VSMALL);
368 
369  forAll(phiCorrIf, facei)
370  {
371  label own = owner[facei];
372  label nei = neighb[facei];
373 
374  psiMaxn[own] = max(psiMaxn[own], psiIf[nei]);
375  psiMinn[own] = min(psiMinn[own], psiIf[nei]);
376 
377  psiMaxn[nei] = max(psiMaxn[nei], psiIf[own]);
378  psiMinn[nei] = min(psiMinn[nei], psiIf[own]);
379 
380  sumPhiBD[own] += phiBDIf[facei];
381  sumPhiBD[nei] -= phiBDIf[facei];
382 
383  scalar phiCorrf = phiCorrIf[facei];
384 
385  if (phiCorrf > 0.0)
386  {
387  sumPhip[own] += phiCorrf;
388  mSumPhim[nei] += phiCorrf;
389  }
390  else
391  {
392  mSumPhim[own] -= phiCorrf;
393  sumPhip[nei] -= phiCorrf;
394  }
395  }
396 
397  forAll(phiCorrBf, patchi)
398  {
399  const fvPatchScalarField& psiPf = psiBf[patchi];
400  const scalarField& phiBDPf = phiBDBf[patchi];
401  const scalarField& phiCorrPf = phiCorrBf[patchi];
402 
403  const labelList& pFaceCells = mesh.boundary()[patchi].faceCells();
404 
405  if (psiPf.coupled())
406  {
407  scalarField psiPNf = psiPf.patchNeighbourField();
408 
409  forAll(phiCorrPf, pFacei)
410  {
411  label pfCelli = pFaceCells[pFacei];
412 
413  psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiPNf[pFacei]);
414  psiMinn[pfCelli] = min(psiMinn[pfCelli], psiPNf[pFacei]);
415  }
416  }
417  else
418  {
419  forAll(phiCorrPf, pFacei)
420  {
421  label pfCelli = pFaceCells[pFacei];
422 
423  psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiPf[pFacei]);
424  psiMinn[pfCelli] = min(psiMinn[pfCelli], psiPf[pFacei]);
425  }
426  }
427 
428  forAll(phiCorrPf, pFacei)
429  {
430  label pfCelli = pFaceCells[pFacei];
431 
432  sumPhiBD[pfCelli] += phiBDPf[pFacei];
433 
434  scalar phiCorrf = phiCorrPf[pFacei];
435 
436  if (phiCorrf > 0.0)
437  {
438  sumPhip[pfCelli] += phiCorrf;
439  }
440  else
441  {
442  mSumPhim[pfCelli] -= phiCorrf;
443  }
444  }
445  }
446 
447  psiMaxn = min(psiMaxn, psiMax);
448  psiMinn = max(psiMinn, psiMin);
449 
450  //scalar smooth = 0.5;
451  //psiMaxn = min((1.0 - smooth)*psiIf + smooth*psiMaxn, psiMax);
452  //psiMinn = max((1.0 - smooth)*psiIf + smooth*psiMinn, psiMin);
453 
454  if (mesh.moving())
455  {
457 
458  psiMaxn =
459  V*((rho/deltaT - Sp)*psiMaxn - Su)
460  - (V0()/deltaT)*rho.oldTime()*psi0
461  + sumPhiBD;
462 
463  psiMinn =
464  V*(Su - (rho/deltaT - Sp)*psiMinn)
465  + (V0/deltaT)*rho.oldTime()*psi0
466  - sumPhiBD;
467  }
468  else
469  {
470  psiMaxn =
471  V*((rho/deltaT - Sp)*psiMaxn - Su - (rho.oldTime()/deltaT)*psi0)
472  + sumPhiBD;
473 
474  psiMinn =
475  V*(Su - (rho/deltaT - Sp)*psiMinn + (rho.oldTime()/deltaT)*psi0)
476  - sumPhiBD;
477  }
478 
479  scalarField sumlPhip(psiIf.size());
480  scalarField mSumlPhim(psiIf.size());
481 
482  for(int j=0; j<nLimiterIter; j++)
483  {
484  sumlPhip = 0.0;
485  mSumlPhim = 0.0;
486 
487  forAll(lambdaIf, facei)
488  {
489  label own = owner[facei];
490  label nei = neighb[facei];
491 
492  scalar lambdaPhiCorrf = lambdaIf[facei]*phiCorrIf[facei];
493 
494  if (lambdaPhiCorrf > 0.0)
495  {
496  sumlPhip[own] += lambdaPhiCorrf;
497  mSumlPhim[nei] += lambdaPhiCorrf;
498  }
499  else
500  {
501  mSumlPhim[own] -= lambdaPhiCorrf;
502  sumlPhip[nei] -= lambdaPhiCorrf;
503  }
504  }
505 
506  forAll(lambdaBf, patchi)
507  {
508  scalarField& lambdaPf = lambdaBf[patchi];
509  const scalarField& phiCorrfPf = phiCorrBf[patchi];
510 
511  const labelList& pFaceCells = mesh.boundary()[patchi].faceCells();
512 
513  forAll(lambdaPf, pFacei)
514  {
515  label pfCelli = pFaceCells[pFacei];
516 
517  scalar lambdaPhiCorrf = lambdaPf[pFacei]*phiCorrfPf[pFacei];
518 
519  if (lambdaPhiCorrf > 0.0)
520  {
521  sumlPhip[pfCelli] += lambdaPhiCorrf;
522  }
523  else
524  {
525  mSumlPhim[pfCelli] -= lambdaPhiCorrf;
526  }
527  }
528  }
529 
530  forAll (sumlPhip, celli)
531  {
532  sumlPhip[celli] =
533  max(min
534  (
535  (sumlPhip[celli] + psiMaxn[celli])/mSumPhim[celli],
536  1.0), 0.0
537  );
538 
539  mSumlPhim[celli] =
540  max(min
541  (
542  (mSumlPhim[celli] + psiMinn[celli])/sumPhip[celli],
543  1.0), 0.0
544  );
545  }
546 
547  const scalarField& lambdam = sumlPhip;
548  const scalarField& lambdap = mSumlPhim;
549 
550  forAll(lambdaIf, facei)
551  {
552  if (phiCorrIf[facei] > 0.0)
553  {
554  lambdaIf[facei] = min
555  (
556  lambdaIf[facei],
557  min(lambdap[owner[facei]], lambdam[neighb[facei]])
558  );
559  }
560  else
561  {
562  lambdaIf[facei] = min
563  (
564  lambdaIf[facei],
565  min(lambdam[owner[facei]], lambdap[neighb[facei]])
566  );
567  }
568  }
569 
570 
571  forAll(lambdaBf, patchi)
572  {
573  fvsPatchScalarField& lambdaPf = lambdaBf[patchi];
574  const scalarField& phiCorrfPf = phiCorrBf[patchi];
575 
576  const labelList& pFaceCells = mesh.boundary()[patchi].faceCells();
577 
578  forAll(lambdaPf, pFacei)
579  {
580  label pfCelli = pFaceCells[pFacei];
581 
582  if (phiCorrfPf[pFacei] > 0.0)
583  {
584  lambdaPf[pFacei] = min(lambdaPf[pFacei], lambdap[pfCelli]);
585  }
586  else
587  {
588  lambdaPf[pFacei] = min(lambdaPf[pFacei], lambdam[pfCelli]);
589  }
590  }
591  }
592 
593  syncTools::syncFaceList(mesh, allLambda, minEqOp<scalar>(), false);
594  }
595 }
596 
597 
598 // ************************ vim: set sw=4 sts=4 et: ************************ //