FreeFOAM The Cross-Platform CFD Toolkit
fvMatrixSolve.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 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
27 
28 template<class Type>
30 (
31  const label patchi,
32  const label facei,
33  const direction cmpt,
34  const scalar value
35 )
36 {
37  if (psi_.needReference())
38  {
39  if (Pstream::master())
40  {
41  internalCoeffs_[patchi][facei].component(cmpt) +=
42  diag()[psi_.mesh().boundary()[patchi].faceCells()[facei]];
43 
44  boundaryCoeffs_[patchi][facei].component(cmpt) +=
45  diag()[psi_.mesh().boundary()[patchi].faceCells()[facei]]
46  *value;
47  }
48  }
49 }
50 
51 
52 template<class Type>
54 (
55  const dictionary& solverControls
56 )
57 {
58  if (debug)
59  {
60  Info<< "fvMatrix<Type>::solve(const dictionary& solverControls) : "
61  "solving fvMatrix<Type>"
62  << endl;
63  }
64 
65  lduMatrix::solverPerformance solverPerfVec
66  (
67  "fvMatrix<Type>::solve",
68  psi_.name()
69  );
70 
71  scalarField saveDiag = diag();
72 
73  Field<Type> source = source_;
74 
75  // At this point include the boundary source from the coupled boundaries.
76  // This is corrected for the implict part by updateMatrixInterfaces within
77  // the component loop.
78  addBoundarySource(source);
79 
80  typename Type::labelType validComponents
81  (
82  pow
83  (
84  psi_.mesh().solutionD(),
86  )
87  );
88 
89  for(direction cmpt=0; cmpt<Type::nComponents; cmpt++)
90  {
91  if (validComponents[cmpt] == -1) continue;
92 
93  // copy field and source
94 
95  scalarField psiCmpt = psi_.internalField().component(cmpt);
96  addBoundaryDiag(diag(), cmpt);
97 
98  scalarField sourceCmpt = source.component(cmpt);
99 
100  FieldField<Field, scalar> bouCoeffsCmpt
101  (
102  boundaryCoeffs_.component(cmpt)
103  );
104 
105  FieldField<Field, scalar> intCoeffsCmpt
106  (
107  internalCoeffs_.component(cmpt)
108  );
109 
110  lduInterfaceFieldPtrsList interfaces =
111  psi_.boundaryField().interfaces();
112 
113  // Use the initMatrixInterfaces and updateMatrixInterfaces to correct
114  // bouCoeffsCmpt for the explicit part of the coupled boundary
115  // conditions
116  initMatrixInterfaces
117  (
118  bouCoeffsCmpt,
119  interfaces,
120  psiCmpt,
121  sourceCmpt,
122  cmpt
123  );
124 
125  updateMatrixInterfaces
126  (
127  bouCoeffsCmpt,
128  interfaces,
129  psiCmpt,
130  sourceCmpt,
131  cmpt
132  );
133 
134  lduMatrix::solverPerformance solverPerf;
135 
136  // Solver call
137  solverPerf = lduMatrix::solver::New
138  (
139  psi_.name() + pTraits<Type>::componentNames[cmpt],
140  *this,
141  bouCoeffsCmpt,
142  intCoeffsCmpt,
143  interfaces,
144  solverControls
145  )->solve(psiCmpt, sourceCmpt, cmpt);
146 
147  solverPerf.print();
148 
149  if
150  (
151  solverPerf.initialResidual() > solverPerfVec.initialResidual()
152  && !solverPerf.singular()
153  )
154  {
155  solverPerfVec = solverPerf;
156  }
157 
158  psi_.internalField().replace(cmpt, psiCmpt);
159  diag() = saveDiag;
160  }
161 
162  psi_.correctBoundaryConditions();
163 
164  return solverPerfVec;
165 }
166 
167 
168 template<class Type>
171 {
172  return solver(psi_.mesh().solverDict(psi_.name()));
173 }
174 
175 template<class Type>
177 {
178  return solve(fvMat_.psi_.mesh().solverDict(fvMat_.psi_.name()));
179 }
180 
181 
182 template<class Type>
184 {
185  return solve(psi_.mesh().solverDict(psi_.name()));
186 }
187 
188 
189 template<class Type>
191 {
192  tmp<Field<Type> > tres(source_);
193  Field<Type>& res = tres();
194 
195  addBoundarySource(res);
196 
197  // Loop over field components
198  for (direction cmpt=0; cmpt<Type::nComponents; cmpt++)
199  {
200  scalarField psiCmpt = psi_.internalField().component(cmpt);
201 
202  scalarField boundaryDiagCmpt(psi_.size(), 0.0);
203  addBoundaryDiag(boundaryDiagCmpt, cmpt);
204 
205  FieldField<Field, scalar> bouCoeffsCmpt
206  (
207  boundaryCoeffs_.component(cmpt)
208  );
209 
210  res.replace
211  (
212  cmpt,
214  (
215  psiCmpt,
216  res.component(cmpt) - boundaryDiagCmpt*psiCmpt,
217  bouCoeffsCmpt,
218  psi_.boundaryField().interfaces(),
219  cmpt
220  )
221  );
222  }
223 
224  return tres;
225 }
226 
227 
228 // ************************ vim: set sw=4 sts=4 et: ************************ //