From 94d3e79a8617f88dc0219cfdeedfa3147833719d Mon Sep 17 00:00:00 2001
From: Adrian Kummerlaender
Date: Mon, 24 Jun 2019 14:43:36 +0200
Subject: Initialize at openlb-1-3
---
src/dynamics/d3q13Helpers.h | 141 ++++++++++++++++++++++++++++++++++++++++++++
1 file changed, 141 insertions(+)
create mode 100644 src/dynamics/d3q13Helpers.h
(limited to 'src/dynamics/d3q13Helpers.h')
diff --git a/src/dynamics/d3q13Helpers.h b/src/dynamics/d3q13Helpers.h
new file mode 100644
index 0000000..24f7a1d
--- /dev/null
+++ b/src/dynamics/d3q13Helpers.h
@@ -0,0 +1,141 @@
+/* This file is part of the OpenLB library
+ *
+ * Copyright (C) 2006, 2007 Jonas Latt
+ * E-mail contact: info@openlb.net
+ * The most recent release of OpenLB can be downloaded at
+ *
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public
+ * License along with this program; if not, write to the Free
+ * Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
+ * Boston, MA 02110-1301, USA.
+*/
+
+/** \file
+ * Helper functions for the implementation of MRT (multiple relaxation
+ * time) dynamics. This file is all about efficiency.
+ * The generic template code is specialized for commonly used Lattices,
+ * so that a maximum performance can be taken out of each case.
+ */
+#ifndef D3Q13_HELPERS_H
+#define D3Q13_HELPERS_H
+
+#include "latticeDescriptors.h"
+#include "core/cell.h"
+#include "core/util.h"
+
+
+namespace olb {
+
+/// Helper functions for the (somewhat special) D3Q13 lattice
+template
+struct d3q13Helpers {
+ typedef descriptors::D3Q13<> DESCRIPTOR;
+
+ /// BGK collision step
+ static T collision( Cell>& cell,
+ T rho, const T u[DESCRIPTOR::d],
+ T lambda_nu, T lambda_nu_prime)
+ {
+ const T lambda_e = descriptors::lambda_e();
+ const T lambda_h = descriptors::lambda_h();
+
+ T uxSqr = util::sqr(u[0]);
+ T uySqr = util::sqr(u[1]);
+ T uzSqr = util::sqr(u[2]);
+ T uSqr = uxSqr + uySqr + uzSqr;
+
+ T s1 = cell[7] + cell[8] + cell[9] + cell[10];
+ T s2 = cell[11] + cell[12];
+ T s3 = cell[1] + cell[2] + cell[3] + cell[4];
+ T s4 = cell[5] + cell[6];
+ T sTot = s1 + s2 + s3 + s4;
+ T d1 = cell[7] + cell[8] - cell[9] - cell[10];
+ T d2 = cell[1] + cell[2] - cell[3] - cell[4];
+
+ T M[13]; // The terms M[0]-M[3] are not used (they correspond
+ // to rho, rho*ux, rho*uy, rho*uz respectively),
+ // but we still use a 13-component vector to preserve
+ // the clarity of the code.
+ M[4] = -(T)12*cell[0] + sTot - (T)11/(T)2;
+ // The 11/2 correction term in M4 accounts for the fact that
+ // cell[i] corresponds to f[i]-ti, and not f[i]
+ M[5] = s1 - (T)2*s2 + s3 - (T)2*s4;
+ M[6] = d1 + d2;
+ M[7] = cell[7] - cell[8] + cell[1] - cell[2];
+ M[8] = cell[11] - cell[12] + cell[5] - cell[6];
+ M[9] = cell[9] - cell[10] + cell[3] - cell[4];
+ M[10] = d1 - d2;
+ M[11] = -cell[7] + cell[8] + s2 + cell[1] - cell[2] - s4;
+ M[12] = cell[9] - cell[10] - cell[11] + cell[12]
+ - cell[3] + cell[4] + cell[5] - cell[6];
+ T Mneq[13]; // The terms Mneq[0]-Mneq[3] are not used (they are
+ // actually all zero, because of conservation laws),
+ // but we still use a 13-component vector to preserve
+ // the clarity of the code.
+ Mneq[4] = M[4] + (T)11/(T)2*rho - (T)13/(T)2*rho*uSqr;
+ Mneq[5] = M[5] - rho*( (T)2*uxSqr-(uySqr+uzSqr) );
+ Mneq[6] = M[6] - rho*( uySqr-uzSqr );
+ Mneq[7] = M[7] - rho*( u[0]*u[1] );
+ Mneq[8] = M[8] - rho*( u[1]*u[2] );
+ Mneq[9] = M[9] - rho*( u[0]*u[2] );
+ Mneq[10] = M[10];
+ Mneq[11] = M[11];
+ Mneq[12] = M[12];
+
+ Mneq[4] *= lambda_e / (T)156;
+ Mneq[5] *= lambda_nu / (T)24;
+ Mneq[6] *= lambda_nu / (T)8;
+ Mneq[7] *= lambda_nu_prime / (T)4;
+ Mneq[8] *= lambda_nu_prime / (T)4;
+ Mneq[9] *= lambda_nu_prime / (T)4;
+ Mneq[10] *= lambda_h / (T)8;
+ Mneq[11] *= lambda_h / (T)8;
+ Mneq[12] *= lambda_h / (T)8;
+
+ T F1 = Mneq[4] + Mneq[5] + Mneq[6];
+ T F2 = Mneq[4] + Mneq[5] - Mneq[6];
+ T F3 = Mneq[4] - (T)2*Mneq[5];
+
+ cell[0] -= (T)-12*Mneq[4];
+ cell[1] -= F1 + Mneq[7] -Mneq[10]+Mneq[11];
+ cell[2] -= F1 - Mneq[7] -Mneq[10]-Mneq[11];
+ cell[3] -= F2 +Mneq[9]+Mneq[10] -Mneq[12];
+ cell[4] -= F2 -Mneq[9]+Mneq[10] +Mneq[12];
+ cell[5] -= F3 +Mneq[8] -Mneq[11]+Mneq[12];
+ cell[6] -= F3 -Mneq[8] -Mneq[11]-Mneq[12];
+ cell[7] -= F1 + Mneq[7] +Mneq[10]-Mneq[11];
+ cell[8] -= F1 - Mneq[7] +Mneq[10]+Mneq[11];
+ cell[9] -= F2 +Mneq[9]-Mneq[10] +Mneq[12];
+ cell[10] -= F2 -Mneq[9]-Mneq[10] -Mneq[12];
+ cell[11] -= F3 +Mneq[8] +Mneq[11]-Mneq[12];
+ cell[12] -= F3 -Mneq[8] +Mneq[11]+Mneq[12];
+
+ return uSqr;
+ }
+
+ /// BGK collision step with density correction
+ static T constRhoCollision( Cell>& cell,
+ T rho, const T u[DESCRIPTOR::d],
+ T ratioRho,
+ T lambda_nu, T lambda_nu_prime)
+ {
+ const T uSqr = util::normSqr(u);
+
+ return uSqr;
+ }
+}; // struct d3q13helpers
+
+} // namespace olb
+
+#endif
--
cgit v1.2.3