OpenLB 1.7
Loading...
Searching...
No Matches
anisoDiscr.h
Go to the documentation of this file.
1/* This file is part of the OpenLB library
2 *
3 * Copyright (C) 2018 Julius Woerner, Albert Mink
4 * E-mail contact: info@openlb.net
5 * The most recent release of OpenLB can be downloaded at
6 * <http://www.openlb.net/>
7 *
8 * This program is free software; you can redistribute it and/or
9 * modify it under the terms of the GNU General Public License
10 * as published by the Free Software Foundation; either version 2
11 * of the License, or (at your option) any later version.
12 *
13 * This program is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public
19 * License along with this program; if not, write to the Free
20 * Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
21 * Boston, MA 02110-1301, USA.
22*/
23
24
31#ifndef ANISO_DISCR_H
32#define ANISO_DISCR_H
33
34#include <iostream>
35#include <vector>
36#include "utilities/omath.h"
37#include <stdlib.h>
38#include <fstream>
39#include <iomanip>
40#include <string>
41#include <math.h>
42#ifdef FEATURE_OPENBLAS
43#include <openblas/lapacke.h>
44#endif
45
47#include "core/util.h"
48
49namespace olb {
50
55double henyeyGreenstein(double cosTheta, double g)
56{
57 return (1-g*g) / util::pow(1+g*g-2*g*cosTheta,1.5);
58}
59
60// check for energy conservation
61template< int q, typename DESCRIPTOR>
62std::array<double,q> testEnergyConservationColumn( const std::array<std::array<double,q>,q>& phi )
63{
64 std::array<double,q> discIntegral;
65 for ( int iVec = 0; iVec < q; iVec++ ) {
66 discIntegral[iVec] = 0;
67 for ( int iSum = 0; iSum < q; iSum++ ) {
68 discIntegral[iVec] += descriptors::t<double,DESCRIPTOR>(iSum) * phi[iSum][iVec];
69 }
70 }
71 return discIntegral;
72}
73
74// check for energy conservation
75template<int q, typename DESCRIPTOR>
76std::array<double,q> testEnergyConservationRow( const std::array<std::array<double,q>,q>& phi )
77{
78 std::array<double,q> discIntegral;
79 for ( int iVec = 0; iVec < q; iVec++ ) {
80 discIntegral[iVec] = 0;
81 for ( int iSum = 0; iSum < q; iSum++ ) {
82 discIntegral[iVec] += descriptors::t<double,DESCRIPTOR>(iSum) * phi[iVec][iSum];
83 }
84 }
85 return discIntegral;
86}
87
88// check for anisotropy conservation
89template< int q>
90std::array<double,q> testAnisotropyConservationColumn( const std::array<std::array<double,q>,q>& phi,
91 const double weights[q], std::array<std::array<double,q>,q>& cosTheta)
92{
93 std::array<double,q> discIntegral;
94 for ( int iVec = 0; iVec < q; iVec++ ) {
95 discIntegral[iVec] = 0;
96 for ( int iSum = 0; iSum < q; iSum++ ) {
97 discIntegral[iVec] += weights[iSum] * cosTheta[iVec][iSum]* phi[iVec][iSum];
98 }
99 }
100 return discIntegral;
101}
102
103
109std::vector<double> linespace( double const stepsize, double const a, double const b )
110{
111 std::vector<double> linspace{}; // initalize to empty
112 if ( util::nearZero( a-b ) ) {
113 return linspace;
114 }
115 int const h = (int) ( (b-a) / stepsize);
116 for (int n = 0; n <= h; n++) {
117 linspace.push_back( a +double(n)/h * b );
118 }
119 return linspace;
120}
121
131template<typename DESCRIPTOR>
132void computeAnisotropyMatrix( double const stepsize, double const anisotropyFactor,
133 double solution[(DESCRIPTOR::q-1)*((DESCRIPTOR::q-1)+1)/2],
134 std::array<std::array<double,DESCRIPTOR::q-1>, DESCRIPTOR::q-1>& phi, int const breakAfter = -1)
135{
136 using namespace descriptors;
137 if ( !DESCRIPTOR::template provides<descriptors::tag::RTLBM>() ) {
138 std::cout << "Warning: Compute anisotropy matrix for wrong latice stencil!" << std::endl;
139 std::cout << "Weight for direction 0 is required to be 0." << std::endl;
140 return;
141 }
142
143 OstreamManager clout( std::cout, "AnisotropyMatrix" );
144 clout << "Call AnistorpiyMatrix()" << std::endl;
145#ifdef FEATURE_OPENBLAS
146 clout << "Compute anisotropy matrix ..." << std::endl;
147 typedef DESCRIPTOR L;
148 int const q = L::q-1;
149 int const mm = 2*q;
150 int const nn = (q+1)*q/2;
151
152 // qxq matrix 0th row and 0th column are empty
153 std::array<std::array<double,q>,q> angleProb;
154 double dotProduct;
155 double normI;
156 double normJ;
157 for (int iPop=0; iPop<q; iPop++) {
158 for (int jPop=0; jPop<q; jPop++) {
159 // shift by 1 due to notation in DESCRIPTOR/DESCRIPTOR
160 // exclude 0th direction
161 dotProduct = c<L>(iPop+1)*c<L>(jPop+1);
162 normI = util::sqrt( c<L>(iPop+1)*c<L>(iPop+1) );
163 normJ = util::sqrt( c<L>(jPop+1)*c<L>(jPop+1) );
164 angleProb[iPop][jPop] = dotProduct / (normI*normJ);
165 }
166 }
167
168 for (int i=0; i<q; i++) {
169 for (int j=0; j<q; j++) {
170 phi[i][j]=1.0;
171 }
172 }
173
174 for ( int i = 0; i < nn; i++ ) {
175 solution[i] = 0;
176 };
177
178 int iter;
179 double U[mm][nn];
180 double U_array[nn*mm];
181 std::vector<double> anisoIterVector = linespace( stepsize, 0, anisotropyFactor );
182
183 // additional condition only for unit testing
184 size_t index = 0;
185 for ( ; index < anisoIterVector.size() && index != (std::size_t)(breakAfter); index++) {
186 // wipe matrices and vectors
187 for (int m = 0; m < mm; m++) {
188 for (int n = 0; n < nn; n++) {
189 U[m][n] = 0;
190 }
191 }
192
193 // compute matrix U, iter handels symmetry
194 iter = 0;
195 for (int m=0; m<q; m++) {
196 for (int n=m; n<q; n++) {
197 U[m][iter] = t<double,L>(m+1)*phi[n][m];
198 U[n][iter] = t<double,L>(n+1)*phi[m][n];
199 U[m+q][iter] = t<double,L>(m+1)*phi[n][m]*angleProb[n][m];
200 U[n+q][iter] = t<double,L>(n+1)*phi[m][n]*angleProb[m][n];
201 iter++;
202 }
203 }
204
205 double sum1;
206 double sum2;
207 // compute vector b
208 for (int n=0; n<q; n++) {
209 // get sum
210 sum1 = 0;
211 sum2 = 0;
212 for (int k=0; k<q; k++) {
213 sum1 += t<double,L>(k+1)*phi[k][n];
214 sum2 += t<double,L>(k+1)*phi[k][n]*angleProb[k][n];
215 }
216 solution[n] = 1 - sum1;
217 solution[q+n] = anisoIterVector[index] - sum2;
218 }
219
220 // transform 2d array to 1d array, column major
221 for ( int n = 0; n < nn; ++n) {
222 for ( int m = 0; m < mm; ++m ) {
223 U_array[n*mm +m] = U[m][n];
224 }
225 }
226
227 //util::print(b, nn, "b");
228
229 // solve Ax = QRx = R'Q'x = b
230 // equiv x = Q*R'\x
231 LAPACKE_dgels( LAPACK_COL_MAJOR, 'N', mm,
232 nn, 1, U_array,
233 mm, solution, nn);
234 //util::print(b, nn, "b after");
235
236 iter = 0;
237 for (int m=0; m<q; m++) {
238 for (int n=m; n<q; n++) {
239 phi[n][m] = phi[n][m]*(1+solution[iter]);
240 phi[m][n] = phi[n][m];
241 iter++;
242 }
243 }
244 //util::print( testEnergyConservationColumn<q>(phi,L::t), "colum sum elementwise" );
245 //util::print( testEnergyConservationRow<q>(phi,L::t), "row sum elementwise" );
246 //util::print( testEnergyConservationSec<q>(phi,L::t,angleProb), "second Moment" );
247 }
248 clout << "Terminate after " << index << " iterations" << std::endl;
249#endif
250}
251
252
253template<typename DESCRIPTOR>
254void computeAnisotropyMatrixKimAndLee( double const anisotropyFactor,
255 std::array<std::array<double,DESCRIPTOR::q>,DESCRIPTOR::q>& phi )
256{
257 OstreamManager clout( std::cout, "AnisotropyMatrix_KimAndLee" );
258 clout << "Compute anisotropy matrix ..." << std::endl;
259 typedef DESCRIPTOR L;
260 int const q = L::q;
261
262 // qxq matrix 0th row and 0th column are empty
263 std::array<std::array<double,q>, q> cosTheta;
264 double dotProduct;
265 double normI;
266 double normJ;
267 for (int iPop=0; iPop<q; iPop++) {
268 for (int jPop=0; jPop<q; jPop++) {
269 dotProduct = descriptors::c<L>(iPop) * descriptors::c<L>(jPop);
270 normI = util::sqrt( util::normSqr<int,3>(descriptors::c<L>(iPop)) );
271 normJ = util::sqrt( util::normSqr<int,3>(descriptors::c<L>(jPop)) );
272 cosTheta[iPop][jPop] = dotProduct /(normI*normJ);
273 if ( util::normSqr<int,3>(descriptors::c<L>(iPop)) == 0 || util::normSqr<int,3>(descriptors::c<L>(jPop)) == 0) {
274 cosTheta[iPop][jPop] = 0.0;
275 }
276 }
277 }
278
279 std::array<std::array<double,q>, q> phaseFunction;
280 for (int m=0; m<q; m++) {
281 for (int n=0; n<q; n++) {
282 phaseFunction[m][n] = henyeyGreenstein(cosTheta[m][n], anisotropyFactor);
283 }
284 }
285
286 for (int m=0; m<q; m++) {
287 for (int n=0; n<q; n++) {
288 dotProduct = 0;
289 for (int i = 0; i < q; ++i) {
290 dotProduct += phaseFunction[m][i]*descriptors::t<double,L>(i);
291 }
292 phi[m][n] = phaseFunction[m][n] / dotProduct;
293 }
294 }
295 //util::print( testEnergyConservationColumn<q>(phi,L::t), "colum sum elementwise" );
296 //util::print( testEnergyConservationRow<q>(phi,L::t), "row sum elementwise" );
297 //util::print( testAnisotropyConservationColumn<q>(phi,L::t,cosTheta), "anisotropyConservation Moment" );
298}
299
300
301} // namespace olb
302
303#endif
class for marking output with some text
cpu::simd::Pack< T > sqrt(cpu::simd::Pack< T > value)
Definition pack.h:100
cpu::simd::Pack< T > pow(cpu::simd::Pack< T > base, cpu::simd::Pack< T > exp)
Definition pack.h:112
bool nearZero(const ADf< T, DIM > &a)
Definition aDiff.h:1087
Top level namespace for all of OpenLB.
std::array< double, q > testEnergyConservationRow(const std::array< std::array< double, q >, q > &phi)
Definition anisoDiscr.h:76
std::vector< double > linespace(double const stepsize, double const a, double const b)
Computes vector [a, a+stepsize, a+2*stepsize, ..., b].
Definition anisoDiscr.h:109
void computeAnisotropyMatrixKimAndLee(double const anisotropyFactor, std::array< std::array< double, DESCRIPTOR::q >, DESCRIPTOR::q > &phi)
Definition anisoDiscr.h:254
double henyeyGreenstein(double cosTheta, double g)
Function to compute Henyey Greenstein phase funtion.
Definition anisoDiscr.h:55
std::array< double, q > testAnisotropyConservationColumn(const std::array< std::array< double, q >, q > &phi, const double weights[q], std::array< std::array< double, q >, q > &cosTheta)
Definition anisoDiscr.h:90
void computeAnisotropyMatrix(double const stepsize, double const anisotropyFactor, double solution[(DESCRIPTOR::q-1) *((DESCRIPTOR::q-1)+1)/2], std::array< std::array< double, DESCRIPTOR::q-1 >, DESCRIPTOR::q-1 > &phi, int const breakAfter=-1)
Definition anisoDiscr.h:132
std::array< double, q > testEnergyConservationColumn(const std::array< std::array< double, q >, q > &phi)
Definition anisoDiscr.h:62
Set of functions commonly used in LB computations – header file.