24#ifndef INAMURO_NEWTON_RAPHSON_DYNAMICS_HH
25#define INAMURO_NEWTON_RAPHSON_DYNAMICS_HH
38template<
typename T,
typename DESCRIPTOR,
typename Dynamics,
typename MOMENTA,
int direction,
int orientation>
40 : legacy::BasicDynamics<T,DESCRIPTOR,MOMENTA>(),
41 _boundaryDynamics(omega),
42 clout(std::cout,
"InamuroNewtonRaphsonDynamics")
44 this->
getName() =
"InamuroNewtonRaphsonDynamics";
46 for (
int iDim = 1; iDim < DESCRIPTOR::d; ++iDim) {
51template<
typename T,
typename DESCRIPTOR,
typename Dynamics,
typename MOMENTA,
int direction,
int orientation>
55 return _boundaryDynamics.computeEquilibrium(iPop, rho, u, uSqr);
58template<
typename T,
typename DESCRIPTOR,
typename Dynamics,
typename MOMENTA,
int direction,
int orientation>
66 MOMENTA().computeRhoU(cell,rho,u);
68 constexpr auto missingIndexes = util::subIndexOutgoing<L,direction,orientation>();
69 std::vector<int> knownIndexes;
71 for (
int iPop = 0; iPop < L::q; ++iPop) {
75 for (
unsigned iPop = 0; iPop < missingIndexes.size(); ++iPop) {
76 test[missingIndexes[iPop]] =
false;
78 for (
int iPop = 0; iPop < L::q; ++iPop) {
80 knownIndexes.push_back(iPop);
84 T approxMomentum[L::d];
86 computeApproxMomentum(approxMomentum,cell,rho,u,_xi,knownIndexes,missingIndexes);
88 T error = computeError(rho, u,approxMomentum);
91 while (error > 1.0e-15) {
94 T gradError[L::d], gradGradError[L::d][L::d];
95 computeGradGradError(gradGradError,gradError,rho,u,_xi,approxMomentum,missingIndexes);
97 bool everythingWentFine = newtonRaphson(_xi,gradError,gradGradError);
98 if ((counter > 100000) || everythingWentFine ==
false) {
103 clout <<
"Failed to converge...." << std::endl;
104 clout <<
"Error = " << error << std::endl;
105 clout <<
"u = (" << rho*u[0];
106 for (
int iD=1; iD<DESCRIPTOR::d; ++iD) {
107 clout <<
", " << rho*u[iD];
109 clout <<
")" << std::endl;
110 clout <<
"uApprox = (" << approxMomentum[0];
111 for (
int iD=1; iD<DESCRIPTOR::d; ++iD) {
112 clout <<
", " << approxMomentum[iD];
114 clout <<
")" << std::endl;
115 clout <<
"xi = (" << _xi[0];
116 for (
int iD=1; iD<DESCRIPTOR::d; ++iD) {
117 clout <<
", " << _xi[iD];
119 clout <<
")" << std::endl;
124 computeApproxMomentum(approxMomentum,cell,rho,u,_xi,knownIndexes,missingIndexes);
125 error = computeError(rho, u,approxMomentum);
131 for (
int iDim = 0; iDim < L::d; ++iDim) {
132 if (direction == iDim) {
137 uCs[iDim] = u[iDim] + _xi[iDim+1-counterDim];
141 T uCsSqr = util::normSqr<T,L::d>(uCs);
143 for (
unsigned iPop = 0; iPop < missingIndexes.size(); ++iPop) {
144 cell[missingIndexes[iPop]] = computeEquilibrium(missingIndexes[iPop],_xi[0],uCs,uCsSqr);
147 _boundaryDynamics.collide(cell, statistics);
150template<
typename T,
typename DESCRIPTOR,
typename Dynamics,
typename MOMENTA,
int direction,
int orientation>
153 return _boundaryDynamics.getOmega();
156template<
typename T,
typename DESCRIPTOR,
typename Dynamics,
typename MOMENTA,
int direction,
int orientation>
159 _boundaryDynamics.setOmega(omega);
162template<
typename T,
typename DESCRIPTOR,
typename Dynamics,
typename MOMENTA,
int direction,
int orientation>
165 const T &rho,
const T u[DESCRIPTOR::d],
const T xi[DESCRIPTOR::d],
166 const std::vector<int> knownIndexes,
const std::vector<int> missingIndexes)
168 typedef DESCRIPTOR L;
172 for (
int iDim = 0; iDim < L::d; ++iDim) {
173 if (direction == iDim) {
175 csVel[iDim] = u[iDim];
178 csVel[iDim] = u[iDim] + xi[iDim+1-counter];
182 T csVelSqr = util::normSqr<T,L::d>(csVel);
184 for (
int iDim = 0; iDim < L::d; ++iDim) {
185 approxMomentum[iDim] = T();
186 for (
unsigned iPop = 0; iPop < knownIndexes.size(); ++iPop) {
187 approxMomentum[iDim] += descriptors::c<L>(knownIndexes[iPop],iDim) *
188 cell[knownIndexes[iPop]];
190 for (
unsigned iPop = 0; iPop < missingIndexes.size(); ++iPop) {
191 approxMomentum[iDim] += descriptors::c<L>(missingIndexes[iPop],iDim) *
192 computeEquilibrium(missingIndexes[iPop],xi[0],csVel,csVelSqr);
197template<
typename T,
typename DESCRIPTOR,
typename Dynamics,
typename MOMENTA,
int direction,
int orientation>
198T InamuroNewtonRaphsonDynamics<T,DESCRIPTOR,Dynamics,MOMENTA,direction,orientation>::
199computeError(
const T &rho,
const T u[DESCRIPTOR::d],
const T approxMomentum[DESCRIPTOR::d])
201 typedef DESCRIPTOR L;
204 for (
int iDim = 0; iDim < L::d; ++iDim) {
205 err += (rho * u[iDim]-approxMomentum[iDim]) * (rho * u[iDim]-approxMomentum[iDim]);
210template<
typename T,
typename DESCRIPTOR,
typename Dynamics,
typename MOMENTA,
int direction,
int orientation>
211void InamuroNewtonRaphsonDynamics<T,DESCRIPTOR,Dynamics,MOMENTA,direction,orientation>::computeGradGradError(
212 T gradGradError[DESCRIPTOR::d][DESCRIPTOR::d],T gradError[DESCRIPTOR::d],
213 const T &rho,
const T u[DESCRIPTOR::d],
const T xi[DESCRIPTOR::d],
214 const T approxMomentum[DESCRIPTOR::d],
215 const std::vector<int> missingIndexes)
217 typedef DESCRIPTOR L;
221 for (
int iDim = 0; iDim < L::d; ++iDim) {
222 if (direction == iDim) {
224 csVel[iDim] = u[iDim];
227 csVel[iDim] = u[iDim] + xi[iDim+1-counter];
230 T csVelSqr = util::normSqr<T,L::d>(csVel);
232 std::vector<T> df[L::d];
234 for (
int iXi = 0; iXi < L::d; ++iXi) {
235 for (
unsigned iPop = 0; iPop < missingIndexes.size(); ++iPop) {
236 df[iXi].push_back(T());
242 for (
unsigned iPop = 0; iPop < missingIndexes.size(); ++iPop) {
244 for (
int iDim = 0; iDim < L::d; ++iDim) {
245 cu += descriptors::c<L>(missingIndexes[iPop],iDim) * csVel[iDim];
247 df[0][iPop] = descriptors::t<T,L>(missingIndexes[iPop])
248 * ((T)1+descriptors::invCs2<T,L>()*cu
249 + 0.5 * descriptors::invCs2<T,L>() * descriptors::invCs2<T,L>() * cu * cu
250 - 0.5 * descriptors::invCs2<T,L>() * csVelSqr);
254 for (
int iDim = 0; iDim < L::d; ++iDim) {
255 if (direction == iDim) {
259 for (
unsigned iPop = 0; iPop < missingIndexes.size(); ++iPop) {
261 for (
int jDim = 0; jDim < L::d; ++jDim) {
262 temp += descriptors::c<L>(missingIndexes[iPop],jDim) * csVel[jDim];
265 df[iDim+1-counter][iPop] = xi[0]*descriptors::t<T,L>(missingIndexes[iPop]) *
266 (descriptors::invCs2<T,L>()*descriptors::c<L>(missingIndexes[iPop],iDim)
267 + descriptors::invCs2<T,L>()*descriptors::invCs2<T,L>()*descriptors::c<L>(missingIndexes[iPop],iDim)*temp
268 - descriptors::invCs2<T,L>()*csVel[iDim]);
273 std::vector<T> ddf[L::d][L::d];
275 for (
int iAlpha = 0; iAlpha < L::d; ++iAlpha) {
276 for (
int iBeta = 0; iBeta < L::d; ++iBeta) {
277 for (
unsigned iPop = 0; iPop < missingIndexes.size(); ++iPop) {
278 ddf[iAlpha][iBeta].push_back(T());
286 for (
int iDim = 0; iDim < L::d; ++iDim) {
287 if (direction == iDim) {
291 for (
unsigned iPop = 0; iPop < missingIndexes.size(); ++iPop) {
293 for (
int jDim = 0; jDim < L::d; ++jDim) {
294 temp += descriptors::c<L>(missingIndexes[iPop],jDim) * csVel[jDim];
297 T d_rho_sa = descriptors::t<T,L>(missingIndexes[iPop]) *
298 (descriptors::invCs2<T,L>()*descriptors::c<L>(missingIndexes[iPop],iDim)
299 + descriptors::invCs2<T,L>()*descriptors::invCs2<T,L>()*descriptors::c<L>(missingIndexes[iPop],iDim)*temp
300 - descriptors::invCs2<T,L>()*csVel[iDim]);
302 ddf[iDim+1-counter][0][iPop] = d_rho_sa;
303 ddf[0][iDim+1-counter][iPop] = d_rho_sa;
308 for (
int iAlpha = 1; iAlpha < L::d; ++iAlpha) {
309 for (
int iBeta = 1; iBeta < L::d; ++iBeta) {
310 for (
unsigned iPop = 0; iPop < missingIndexes.size(); ++iPop) {
311 ddf[iAlpha][iBeta][iPop] = descriptors::t<T,L>(missingIndexes[iPop])*xi[0] *
312 descriptors::invCs2<T,L>()*descriptors::invCs2<T,L>()*descriptors::c<L>(missingIndexes[iPop],iAlpha)*descriptors::c<L>(missingIndexes[iPop],iBeta);
313 if (iAlpha == iBeta) {
314 ddf[iAlpha][iBeta][iPop] -= descriptors::t<T,L>(missingIndexes[iPop])*xi[0] * descriptors::invCs2<T,L>();
322 for (
int iDim = 0; iDim < L::d; ++iDim) {
324 gradError[iDim] = T();
325 for (
int jDim = 0; jDim < L::d; ++jDim) {
327 for (
unsigned iPop = 0; iPop < missingIndexes.size(); ++iPop) {
328 du[jDim] += descriptors::c<L>(missingIndexes[iPop],jDim) * df[iDim][iPop];
330 gradError[iDim] += (approxMomentum[jDim]- rho * u[jDim]) * du[jDim];
332 gradError[iDim] *= (T)2;
337 for (
int iAlpha = 0; iAlpha < L::d; ++iAlpha) {
338 for (
int iBeta = 0; iBeta < L::d; ++iBeta) {
339 gradGradError[iAlpha][iBeta] = T();
341 T duAlpha[L::d], duBeta[L::d], ddu[L::d];
342 for (
int iDim = 0; iDim < L::d; ++iDim) {
346 for (
unsigned iPop = 0; iPop < missingIndexes.size(); ++iPop) {
347 duAlpha[iDim] += descriptors::c<L>(missingIndexes[iPop],iDim)
350 duBeta[iDim] += descriptors::c<L>(missingIndexes[iPop],iDim)
353 ddu[iDim] += descriptors::c<L>(missingIndexes[iPop],iDim)
354 * ddf[iAlpha][iBeta][iPop];
356 gradGradError[iAlpha][iBeta] += duAlpha[iDim] * duBeta[iDim]
357 + (approxMomentum[iDim]-rho * u[iDim]) * ddu[iDim];
359 gradGradError[iAlpha][iBeta] *= (T)2;
366template<
typename T,
typename DESCRIPTOR,
typename Dynamics,
typename MOMENTA,
int direction,
int orientation>
367bool InamuroNewtonRaphsonDynamics<T,DESCRIPTOR,Dynamics,MOMENTA,direction,orientation>::
368newtonRaphson(T xi[DESCRIPTOR::d],
369 const T gradError[DESCRIPTOR::d],
370 const T gradGradError[DESCRIPTOR::d][DESCRIPTOR::d])
372 typedef DESCRIPTOR L;
374 T invGradGradError[L::d][L::d];
375 bool inversion = invert(gradGradError,invGradGradError);
377 for (
int iXi = 0; iXi < L::d; ++iXi) {
379 for (
int iAlpha = 0; iAlpha < L::d; ++iAlpha) {
380 correction += invGradGradError[iXi][iAlpha] * gradError[iAlpha];
383 xi[iXi] -= correction;
389template<
typename T,
typename DESCRIPTOR,
typename Dynamics,
typename MOMENTA,
int direction,
int orientation>
390bool InamuroNewtonRaphsonDynamics<T,DESCRIPTOR,Dynamics,MOMENTA,direction,orientation>::
391invert(
const T a[2][2],T invA[2][2])
393 T detA = a[0][0]*a[1][1] - a[0][1]*a[1][0];
396 clout <<
"error detA too small! = " << detA << std::endl;
397 for (
int iAlpha = 0; iAlpha < 2; ++iAlpha) {
398 for (
int iBeta = 0; iBeta < 2; ++iBeta) {
399 clout << a[iAlpha][iBeta] <<
"\t";
406 invA[0][0] = a[1][1];
407 invA[1][1] = a[0][0];
409 invA[0][1] = -a[1][0];
410 invA[1][0] = -a[0][1];
412 for (
int iPop = 0; iPop < 2; ++iPop) {
413 for (
int jPop = 0; jPop < 2; ++jPop) {
414 invA[iPop][jPop] /= detA;
421template<
typename T,
typename DESCRIPTOR,
typename Dynamics,
typename MOMENTA,
int direction,
int orientation>
422bool InamuroNewtonRaphsonDynamics<T,DESCRIPTOR,Dynamics,MOMENTA,direction,orientation>::invert(
const T a[3][3],T invA[3][3])
424 T detA = a[0][0]*a[1][1]*a[2][2] + a[1][0]*a[2][1]*a[0][2] + a[2][0]*a[0][1]*a[1][2]
425 - a[0][0]*a[2][1]*a[1][2] - a[2][0]*a[1][1]*a[0][2] - a[1][0]*a[0][1]*a[2][2];
429 clout <<
"Error: detA too small! = " << detA << std::endl;
430 for (
int iAlpha = 0; iAlpha < 3; ++iAlpha) {
431 for (
int iBeta = 0; iBeta < 3; ++iBeta) {
432 clout << a[iAlpha][iBeta] <<
"\t";
439 invA[0][0] = a[1][1]*a[2][2]-a[1][2]*a[2][1];
440 invA[0][1] = a[0][2]*a[2][1]-a[0][1]*a[2][2];
441 invA[0][2] = a[0][1]*a[1][2]-a[0][2]*a[1][1];
443 invA[1][0] = a[1][2]*a[2][0]-a[1][0]*a[2][2];
444 invA[1][1] = a[0][0]*a[2][2]-a[0][2]*a[2][0];
445 invA[1][2] = a[0][2]*a[1][0]-a[0][0]*a[1][2];
447 invA[2][0] = a[1][0]*a[2][1]-a[1][1]*a[2][0];
448 invA[2][1] = a[0][1]*a[2][0]-a[0][0]*a[2][1];
449 invA[2][2] = a[0][0]*a[1][1]-a[0][1]*a[1][0];
451 for (
int iPop = 0; iPop < 3; ++iPop) {
452 for (
int jPop = 0; jPop < 3; ++jPop) {
453 invA[iPop][jPop] /= detA;
Highest-level interface to Cell data.
Highest-level interface to read-only Cell data.
This class computes the inamuro BC with general dynamics.
T computeEquilibrium(int iPop, T rho, const T u[DESCRIPTOR::d], T uSqr) const override
Compute equilibrium distribution function.
T getOmega() const
Get local relaxation parameter of the dynamics.
void setOmega(T omega)
Set local relaxation parameter of the dynamics.
CellStatistic< T > collide(Cell< T, DESCRIPTOR > &cell, LatticeStatistics< T > &statistics) override
Collision step.
InamuroNewtonRaphsonDynamics(T omega)
Constructor.
Descriptor for all types of 2D and 3D lattices.
cpu::simd::Pack< T > sqrt(cpu::simd::Pack< T > value)
cpu::simd::Pack< T > fabs(cpu::simd::Pack< T > value)
Top level namespace for all of OpenLB.
Return value of any collision.
virtual std::string getName() const
Return human-readable name.
Set of functions commonly used in LB computations – header file.