Feature: block coupling rewrite

This commit is contained in:
Hrvoje Jasak 2015-04-27 18:12:34 +01:00 committed by Dominik Christ
parent 2e0286b2c1
commit 75259a047b
2 changed files with 14 additions and 54 deletions

View file

@ -38,13 +38,8 @@ Description
#include "fieldTypes.H"
#include "Time.H"
#include "fvMesh.H"
#include "fvBlockMatrix.H"
#include "blockLduSolvers.H"
#include "VectorNFieldTypes.H"
#include "volVectorNFields.H"
#include "blockVectorNMatrices.H"
#include "blockMatrixTools.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -92,53 +87,27 @@ int main(int argc, char *argv[])
TsEqn.relax();
// Prepare block system
BlockLduMatrix<vector2> blockM(mesh);
//- Transfer the coupled interface list for processor/cyclic/etc.
// boundaries
blockM.interfaces() = blockT.boundaryField().blockInterfaces();
// Grab block diagonal and set it to zero
Field<tensor2>& d = blockM.diag().asSquare();
d = tensor2::zero;
// Grab linear off-diagonal and set it to zero
Field<vector2>& l = blockM.lower().asLinear();
Field<vector2>& u = blockM.upper().asLinear();
u = vector2::zero;
l = vector2::zero;
vector2Field& blockX = blockT.internalField();
vector2Field blockB(mesh.nCells(), vector2::zero);
fvBlockMatrix<vector2> blockM(blockT);
//- Inset equations into block Matrix
blockMatrixTools::insertEquation(0, TEqn, blockM, blockX, blockB);
blockMatrixTools::insertEquation(1, TsEqn, blockM, blockX, blockB);
blockM.insertEquation(0, TEqn);
blockM.insertEquation(1, TsEqn);
//- Add off-diagonal terms and remove from block source
forAll(d, i)
{
d[i](0, 1) = -alpha.value()*mesh.V()[i];
d[i](1, 0) = -alpha.value()*mesh.V()[i];
//- Add off-diagonal coupling terms
scalarField coupling(mesh.nCells(), -alpha.value());
blockB[i][0] -= alpha.value()*blockX[i][1]*mesh.V()[i];
blockB[i][1] -= alpha.value()*blockX[i][0]*mesh.V()[i];
}
blockM.insertEquationCoupling(0, 1, coupling);
blockM.insertEquationCoupling(1, 0, coupling);
// Update source coupling: coupling terms eliminated from source
blockM.updateSourceCoupling();
//- Block coupled solver call
BlockSolverPerformance<vector2> solverPerf =
BlockLduSolver<vector2>::New
(
blockT.name(),
blockM,
mesh.solutionDict().solver(blockT.name())
)->solve(blockX, blockB);
solverPerf.print();
blockM.solve();
// Retrieve solution
blockMatrixTools::blockRetrieve(0, T.internalField(), blockX);
blockMatrixTools::blockRetrieve(1, Ts.internalField(), blockX);
blockM.retrieveSolution(0, T.internalField());
blockM.retrieveSolution(1, Ts.internalField());
T.correctBoundaryConditions();
Ts.correctBoundaryConditions();

View file

@ -44,15 +44,6 @@
dimensionedVector2("zero", dimless, vector2::zero)
);
{
vector2Field& blockX = blockT.internalField();
blockMatrixTools::blockInsert(0, T.internalField(), blockX);
blockMatrixTools::blockInsert(1, Ts.internalField(), blockX);
}
blockT.write();
Info<< "Reading field U\n" << endl;
volVectorField U