Skip to content

Commit 0b9018f

Browse files
LeanFangLeanFang
authored
Irk pimple (#782)
* IrkPimple * format * relaxStage = 0.8 * IRKDict optional * added IRK pimple reg test (laminar) * reg test UNorm ref updated * unitTest
1 parent 3a6d290 commit 0b9018f

9 files changed

+538
-0
lines changed

‎src/adjoint/DASolver/DAIrkPimpleFoam/DAIrkPimpleFoam.C

+265
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,265 @@
1+
/*---------------------------------------------------------------------------*\
2+
3+
DAFoam : Discrete Adjoint with OpenFOAM
4+
Version : v4
5+
6+
This class is modified from OpenFOAM's source code
7+
applications/solvers/incompressible/pimpleFoam
8+
9+
OpenFOAM: The Open Source CFD Toolbox
10+
11+
Copyright (C): 2011-2016 OpenFOAM Foundation
12+
13+
OpenFOAM License:
14+
15+
OpenFOAM is free software: you can redistribute it and/or modify it
16+
under the terms of the GNU General Public License as published by
17+
the Free Software Foundation, either version 3 of the License, or
18+
(at your option) any later version.
19+
20+
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
21+
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
22+
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
23+
for more details.
24+
25+
You should have received a copy of the GNU General Public License
26+
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
27+
28+
\*---------------------------------------------------------------------------*/
29+
30+
#include "DAIrkPimpleFoam.H"
31+
32+
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
33+
34+
namespace Foam
35+
{
36+
37+
defineTypeNameAndDebug(DAIrkPimpleFoam, 0);
38+
addToRunTimeSelectionTable(DASolver, DAIrkPimpleFoam, dictionary);
39+
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
40+
41+
DAIrkPimpleFoam::DAIrkPimpleFoam(
42+
char* argsAll,
43+
PyObject* pyOptions)
44+
: DASolver(argsAll, pyOptions)
45+
{
46+
}
47+
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
48+
49+
void DAIrkPimpleFoam::initSolver()
50+
{
51+
/*
52+
Description:
53+
Initialize variables for DASolver
54+
*/
55+
daOptionPtr_.reset(new DAOption(meshPtr_(), pyOptions_));
56+
}
57+
58+
label DAIrkPimpleFoam::solvePrimal()
59+
{
60+
/*
61+
Description:
62+
Call the primal solver to get converged state variables
63+
64+
Output:
65+
state variable vector
66+
*/
67+
68+
Foam::argList& args = argsPtr_();
69+
#include "createTime.H"
70+
#include "createMesh.H"
71+
#include "initContinuityErrs.H"
72+
#include "createControl.H"
73+
#include "createFieldsIrkPimple.H"
74+
#include "CourantNo.H"
75+
76+
// Turbulence disabled
77+
//turbulence->validate();
78+
79+
// Get nu from memory
80+
volScalarField& nu = const_cast<volScalarField&>(mesh.thisDb().lookupObject<volScalarField>("nu"));
81+
82+
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
83+
84+
Info << "\nStarting time loop\n"
85+
<< endl;
86+
87+
// Radau23 coefficients and weights
88+
scalar D10 = -2;
89+
scalar D11 = 3.0 / 2;
90+
scalar D12 = 1.0 / 2;
91+
scalar D20 = 2;
92+
scalar D21 = -9.0 / 2;
93+
scalar D22 = 5.0 / 2;
94+
scalar w1 = 3.0 / 4;
95+
scalar w2 = 1.0 / 4;
96+
97+
// get IRKDict settings, default to Radau23 for now
98+
IOdictionary IRKDict(
99+
IOobject(
100+
"IRKDict",
101+
mesh.time().system(),
102+
mesh,
103+
IOobject::READ_IF_PRESENT,
104+
IOobject::NO_WRITE));
105+
106+
scalar relaxU = 1.0;
107+
if (IRKDict.found("relaxU"))
108+
{
109+
if (IRKDict.getScalar("relaxU") > 0)
110+
{
111+
relaxU = IRKDict.getScalar("relaxU");
112+
}
113+
}
114+
115+
scalar relaxP = 1.0;
116+
if (IRKDict.found("relaxP"))
117+
{
118+
if (IRKDict.getScalar("relaxP") > 0)
119+
{
120+
relaxP = IRKDict.getScalar("relaxP");
121+
}
122+
}
123+
124+
scalar relaxPhi = 1.0;
125+
if (IRKDict.found("relaxPhi"))
126+
{
127+
if (IRKDict.getScalar("relaxPhi") > 0)
128+
{
129+
relaxPhi = IRKDict.getScalar("relaxPhi");
130+
}
131+
}
132+
133+
scalar relaxStage1 = 0.8;
134+
if (IRKDict.found("relaxStage1"))
135+
{
136+
if (IRKDict.getScalar("relaxStage1") > 0)
137+
{
138+
relaxStage1 = IRKDict.getScalar("relaxStage1");
139+
}
140+
}
141+
142+
scalar relaxStage2 = 0.8;
143+
if (IRKDict.found("relaxStage2"))
144+
{
145+
if (IRKDict.getScalar("relaxStage2") > 0)
146+
{
147+
relaxStage2 = IRKDict.getScalar("relaxStage2");
148+
}
149+
}
150+
151+
scalar relaxUEqn = 1.0;
152+
153+
label maxSweep = 10;
154+
if (IRKDict.found("maxSweep"))
155+
{
156+
if (IRKDict.getLabel("maxSweep") > 0)
157+
{
158+
maxSweep = IRKDict.getLabel("maxSweep");
159+
}
160+
}
161+
162+
// Duplicate state variables for stages
163+
volVectorField U1("U1", U);
164+
volVectorField U2("U2", U);
165+
volScalarField p1("p1", p);
166+
volScalarField p2("p2", p);
167+
surfaceScalarField phi1("phi1", phi);
168+
surfaceScalarField phi2("phi2", phi);
169+
170+
// Settings for stage pressure
171+
mesh.setFluxRequired(p1.name());
172+
mesh.setFluxRequired(p2.name());
173+
174+
// Note: below is not working somehow...
175+
/*
176+
// IO settings for internal stages
177+
U1.writeOpt() = IOobject::AUTO_WRITE;
178+
p1.writeOpt() = IOobject::AUTO_WRITE;
179+
phi1.writeOpt() = IOobject::AUTO_WRITE;
180+
*/
181+
182+
// Initialize oldTime() for under-relaxation
183+
volVectorField U1OldTime("U1OldTime", U1);
184+
volVectorField U2OldTime("U2OldTime", U2);
185+
volScalarField p1OldTime("p1OldTime", p1);
186+
volScalarField p2OldTime("p2OldTime", p2);
187+
surfaceScalarField phi1OldTime("phi1OldTime", phi1);
188+
surfaceScalarField phi2OldTime("phi2OldTime", phi2);
189+
190+
// Numerical settings
191+
word divUScheme = "div(phi,U)";
192+
//word divGradUScheme = "div((nuEff*dev2(T(grad(U)))))";
193+
194+
const fvSolution& myFvSolution = mesh.thisDb().lookupObject<fvSolution>("fvSolution");
195+
dictionary solverDictU = myFvSolution.subDict("solvers").subDict("U");
196+
dictionary solverDictP = myFvSolution.subDict("solvers").subDict("p");
197+
198+
while (runTime.run())
199+
{
200+
201+
#include "CourantNo.H"
202+
203+
++runTime;
204+
205+
Info << "Time = " << runTime.timeName() << nl << endl;
206+
207+
scalar deltaT = runTime.deltaTValue();
208+
209+
// --- GS sweeps for IRK-PIMPLE
210+
label sweepIndex = 0;
211+
while (sweepIndex < maxSweep)
212+
{
213+
Info << "Block GS sweep = " << sweepIndex + 1 << endl;
214+
215+
{
216+
#include "U1EqnIrkPimple.H"
217+
218+
while (pimple.correct())
219+
{
220+
#include "p1EqnIrkPimple.H"
221+
}
222+
}
223+
224+
{
225+
#include "U2EqnIrkPimple.H"
226+
227+
while (pimple.correct())
228+
{
229+
#include "p2EqnIrkPimple.H"
230+
}
231+
}
232+
233+
sweepIndex++;
234+
}
235+
236+
// Update new step values before write-to-disk
237+
U = U2;
238+
p = p2;
239+
phi = phi2;
240+
241+
runTime.write();
242+
// Also write internal stages to disk (Radau23)
243+
U1.write();
244+
p1.write();
245+
phi1.write();
246+
247+
// Use old step as initial guess for the next step
248+
U1 = U;
249+
U1.correctBoundaryConditions();
250+
p1 = p;
251+
p1.correctBoundaryConditions();
252+
phi1 = phi;
253+
254+
runTime.printExecutionTime(Info);
255+
}
256+
257+
Info << "End\n"
258+
<< endl;
259+
260+
return 0;
261+
}
262+
263+
} // End namespace Foam
264+
265+
// ************************************************************************* //

‎src/adjoint/DASolver/DAIrkPimpleFoam/DAIrkPimpleFoam.H

+86
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,86 @@
1+
/*---------------------------------------------------------------------------*\
2+
3+
DAFoam : Discrete Adjoint with OpenFOAM
4+
Version : v4
5+
6+
Description:
7+
Child class for DAIrkPimpleFoam
8+
9+
This class is modified from OpenFOAM's source code
10+
applications/solvers/incompressible/pimpleFoam
11+
12+
OpenFOAM: The Open Source CFD Toolbox
13+
14+
Copyright (C): 2011-2016 OpenFOAM Foundation
15+
16+
OpenFOAM License:
17+
18+
OpenFOAM is free software: you can redistribute it and/or modify it
19+
under the terms of the GNU General Public License as published by
20+
the Free Software Foundation, either version 3 of the License, or
21+
(at your option) any later version.
22+
23+
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
24+
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
25+
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
26+
for more details.
27+
28+
You should have received a copy of the GNU General Public License
29+
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
30+
31+
\*---------------------------------------------------------------------------*/
32+
33+
#ifndef DAIrkPimpleFoam_H
34+
#define DAIrkPimpleFoam_H
35+
36+
#include "DASolver.H"
37+
#include "addToRunTimeSelectionTable.H"
38+
#include "singlePhaseTransportModel.H"
39+
#include "turbulentTransportModel.H"
40+
//#include "pimpleControlDF.H" // A modified version of pimpleControl, Basically, we disable the output to the screen
41+
#include "pimpleControl.H" // origional pimpleControl
42+
43+
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
44+
45+
namespace Foam
46+
{
47+
48+
/*---------------------------------------------------------------------------*\
49+
Class DAIrkPimpleFoam Declaration
50+
\*---------------------------------------------------------------------------*/
51+
52+
class DAIrkPimpleFoam
53+
: public DASolver
54+
{
55+
56+
protected:
57+
public:
58+
TypeName("DAIrkPimpleFoam");
59+
// Constructors
60+
61+
//- Construct from components
62+
DAIrkPimpleFoam(
63+
char* argsAll,
64+
PyObject* pyOptions);
65+
66+
//- Destructor
67+
virtual ~DAIrkPimpleFoam()
68+
{
69+
}
70+
71+
/// initialize fields and variables
72+
virtual void initSolver();
73+
74+
/// solve the primal equations
75+
virtual label solvePrimal();
76+
};
77+
78+
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
79+
80+
} // End namespace Foam
81+
82+
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
83+
84+
#endif
85+
86+
// ************************************************************************* //

‎src/adjoint/DASolver/DAIrkPimpleFoam/U1EqnIrkPimple.H

+24
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,24 @@
1+
// Initialize U1Eqn w/o ddt term
2+
fvVectorMatrix U1Eqn(
3+
fvm::div(phi1, U1, divUScheme)
4+
//+ turbulence->divDevReff(U)
5+
- fvm::laplacian(nu, U1));
6+
7+
// Update U1Eqn with pseudo-spectral terms
8+
forAll(U1, cellI)
9+
{
10+
scalar meshV = U1.mesh().V()[cellI];
11+
12+
// Add D11 / halfDeltaT[i] * V() to diagonal
13+
U1Eqn.diag()[cellI] += D11 / deltaT * meshV; // Use one seg for now: halfDeltaTList[segI]
14+
15+
// Minus D10 / halfDeltaT[i] * T0 * V() to source term
16+
U1Eqn.source()[cellI] -= D10 / deltaT * U[cellI] * meshV;
17+
18+
// Minus D12 / halfDeltaT[i] * T2 * V() to source term
19+
U1Eqn.source()[cellI] -= D12 / deltaT * U2[cellI] * meshV;
20+
}
21+
22+
U1Eqn.relax(relaxUEqn);
23+
24+
solve(U1Eqn == -fvc::grad(p1), solverDictU);

0 commit comments

Comments
 (0)