OpenLB 1.7
Loading...
Searching...
No Matches
freeEnergyPostProcessor3D.hh
Go to the documentation of this file.
1/* This file is part of the OpenLB library
2 *
3 * Copyright (C) 2018 Sam Avis, Robin Trunk
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#ifndef FREE_ENERGY_POST_PROCESSOR_3D_HH
25#define FREE_ENERGY_POST_PROCESSOR_3D_HH
26
27#include "freeEnergyPostProcessor3D.h"
28
29namespace olb {
30
32
33template<typename T, typename DESCRIPTOR>
34FreeEnergyChemicalPotentialCoupling3D <T,DESCRIPTOR>::FreeEnergyChemicalPotentialCoupling3D (
35 int x0_, int x1_, int y0_, int y1_, int z0_, int z1_, T alpha_, T kappa1_, T kappa2_, T kappa3_,
36 std::vector<BlockStructureD<3>*> partners_)
37 : x0(x0_), x1(x1_), y0(y0_), y1(y1_), z0(z0_), z1(z1_), alpha(alpha_),
38 kappa1(kappa1_), kappa2(kappa2_), kappa3(kappa3_), partners(partners_)
39{
40 this->getName() = "FreeEnergyChemicalPotentialCoupling3D";
41}
42
43template<typename T, typename DESCRIPTOR>
44FreeEnergyChemicalPotentialCoupling3D <T,DESCRIPTOR>::FreeEnergyChemicalPotentialCoupling3D (
45 T alpha_, T kappa1_, T kappa2_, T kappa3_,
46 std::vector<BlockStructureD<3>*> partners_)
47 : x0(0), x1(0), y0(0), y1(0), z0(0), z1(0), alpha(alpha_),
48 kappa1(kappa1_), kappa2(kappa2_), kappa3(kappa3_), partners(partners_)
49{
50 this->getName() = "FreeEnergyChemicalPotentialCoupling3D";
51}
52
53template<typename T, typename DESCRIPTOR>
54void FreeEnergyChemicalPotentialCoupling3D<T,DESCRIPTOR>::processSubDomain (
55 BlockLattice<T,DESCRIPTOR>& blockLattice,
56 int x0_, int x1_, int y0_, int y1_, int z0_, int z1_ )
57{
58 // If partners.size() == 1: two fluid components
59 // If partners.size() == 2: three fluid components
60 BlockLattice<T,DESCRIPTOR> *partnerLattice1 = static_cast<BlockLattice<T,DESCRIPTOR> *>(partners[0]);
61 BlockLattice<T,DESCRIPTOR> *partnerLattice2 = 0;
62 if (partners.size() > 1) {
63 partnerLattice2 = static_cast<BlockLattice<T,DESCRIPTOR> *>(partners[1]);
64 }
65
66 int newX0, newX1, newY0, newY1, newZ0, newZ1;
67 if ( util::intersect ( x0, x1, y0, y1, z0, z1,
68 x0_, x1_, y0_, y1_, z0, z1,
69 newX0, newX1, newY0, newY1, newZ0, newZ1 ) ) {
70
71 // compute the density fields for each lattice
72 auto& rhoField = blockLattice.template getField<RHO_CACHE>();
73
74 for (int iX=newX0-1; iX<=newX1+1; ++iX)
75 for (int iY=newY0-1; iY<=newY1+1; ++iY)
76 for (int iZ=newZ0-1; iZ<=newZ1+1; ++iZ) {
77 rhoField[0][blockLattice.getCellId(iX, iY, iZ)] = blockLattice.get(iX,iY,iZ).computeRho();
78 }
79 for (int iX=newX0-1; iX<=newX1+1; ++iX)
80 for (int iY=newY0-1; iY<=newY1+1; ++iY)
81 for (int iZ=newZ0-1; iZ<=newZ1+1; ++iZ) {
82 rhoField[1][blockLattice.getCellId(iX, iY, iZ)] = partnerLattice1->get(iX,iY,iZ).computeRho();
83 }
84 if (partners.size() > 1) {
85 for (int iX=newX0-1; iX<=newX1+1; ++iX)
86 for (int iY=newY0-1; iY<=newY1+1; ++iY)
87 for (int iZ=newZ0-1; iZ<=newZ1+1; ++iZ) {
88 rhoField[2][blockLattice.getCellId(iX, iY, iZ)] = partnerLattice2->get(iX,iY,iZ).computeRho();
89 }
90 }
91
92 // calculate chemical potential
93 for (int iX=newX0; iX<=newX1; ++iX) {
94 for (int iY=newY0; iY<=newY1; ++iY) {
95 for (int iZ=newZ0; iZ<=newZ1; ++iZ) {
96 T densitySum = rhoField[0][blockLattice.getCellId(iX, iY, iZ)]
97 + rhoField[1][blockLattice.getCellId(iX, iY, iZ)];
98 T densityDifference = rhoField[0][blockLattice.getCellId(iX, iY, iZ)]
99 - rhoField[1][blockLattice.getCellId(iX, iY, iZ)];
100 if (partners.size() > 1) {
101 densitySum -= rhoField[2][blockLattice.getCellId(iX, iY, iZ)];
102 densityDifference -= rhoField[2][blockLattice.getCellId(iX, iY, iZ)];
103 }
104 T term1 = 0.125 * kappa1 * (densitySum)
105 * (densitySum-1.) * (densitySum-2.);
106 T term2 = 0.125 * kappa2 * (densityDifference)
107 * (densityDifference-1.) * (densityDifference-2.);
108 T term3 = 0.;
109 if (partners.size() > 1) {
110 T rho3 = rhoField[2][blockLattice.getCellId(iX, iY, iZ)];
111 term3 = kappa3 * rho3 * (rho3 - 1.) * (2.*rho3 - 1.);
112 }
113
114 T laplaceRho1 = 1.0 / 6.0 * (
115 rhoField[0][blockLattice.getCellId(iX, iY-1, iZ-1)]
116 + rhoField[0][blockLattice.getCellId(iX-1, iY, iZ-1)]
117 + 2. * rhoField[0][blockLattice.getCellId(iX, iY, iZ-1)]
118 + rhoField[0][blockLattice.getCellId(iX+1, iY, iZ-1)]
119 + rhoField[0][blockLattice.getCellId(iX, iY+1, iZ-1)]
120 + rhoField[0][blockLattice.getCellId(iX-1, iY-1, iZ)]
121 + 2. * rhoField[0][blockLattice.getCellId(iX, iY-1, iZ)]
122 + rhoField[0][blockLattice.getCellId(iX+1, iY-1, iZ)]
123 + 2. * rhoField[0][blockLattice.getCellId(iX-1, iY, iZ)]
124 -24. * rhoField[0][blockLattice.getCellId(iX, iY, iZ)]
125 + 2. * rhoField[0][blockLattice.getCellId(iX+1, iY, iZ)]
126 + rhoField[0][blockLattice.getCellId(iX-1, iY+1, iZ)]
127 + 2. * rhoField[0][blockLattice.getCellId(iX, iY+1, iZ)]
128 + rhoField[0][blockLattice.getCellId(iX+1, iY+1, iZ)]
129 + rhoField[0][blockLattice.getCellId(iX, iY-1, iZ+1)]
130 + rhoField[0][blockLattice.getCellId(iX-1, iY, iZ+1)]
131 + 2. * rhoField[0][blockLattice.getCellId(iX, iY, iZ+1)]
132 + rhoField[0][blockLattice.getCellId(iX+1, iY, iZ+1)]
133 + rhoField[0][blockLattice.getCellId(iX, iY+1, iZ+1)]
134 );
135
136 T laplaceRho2 = 1.0 / 6.0 * (
137 rhoField[1][blockLattice.getCellId(iX, iY-1, iZ-1)]
138 + rhoField[1][blockLattice.getCellId(iX-1, iY, iZ-1)]
139 + 2. * rhoField[1][blockLattice.getCellId(iX, iY, iZ-1)]
140 + rhoField[1][blockLattice.getCellId(iX+1, iY, iZ-1)]
141 + rhoField[1][blockLattice.getCellId(iX, iY+1, iZ-1)]
142 + rhoField[1][blockLattice.getCellId(iX-1, iY-1, iZ)]
143 + 2. * rhoField[1][blockLattice.getCellId(iX, iY-1, iZ)]
144 + rhoField[1][blockLattice.getCellId(iX+1, iY-1, iZ)]
145 + 2. * rhoField[1][blockLattice.getCellId(iX-1, iY, iZ)]
146 -24. * rhoField[1][blockLattice.getCellId(iX, iY, iZ)]
147 + 2. * rhoField[1][blockLattice.getCellId(iX+1, iY, iZ)]
148 + rhoField[1][blockLattice.getCellId(iX-1, iY+1, iZ)]
149 + 2. * rhoField[1][blockLattice.getCellId(iX, iY+1, iZ)]
150 + rhoField[1][blockLattice.getCellId(iX+1, iY+1, iZ)]
151 + rhoField[1][blockLattice.getCellId(iX, iY-1, iZ+1)]
152 + rhoField[1][blockLattice.getCellId(iX-1, iY, iZ+1)]
153 + 2. * rhoField[1][blockLattice.getCellId(iX, iY, iZ+1)]
154 + rhoField[1][blockLattice.getCellId(iX+1, iY, iZ+1)]
155 + rhoField[1][blockLattice.getCellId(iX, iY+1, iZ+1)]
156 );
157
158 T laplaceRho3 = 0.;
159 if (partners.size() > 1) {
160 laplaceRho3 = 1.0 / 6.0 * (
161 rhoField[2][blockLattice.getCellId(iX, iY-1, iZ-1)]
162 + rhoField[2][blockLattice.getCellId(iX-1, iY, iZ-1)]
163 + 2. * rhoField[2][blockLattice.getCellId(iX, iY, iZ-1)]
164 + rhoField[2][blockLattice.getCellId(iX+1, iY, iZ-1)]
165 + rhoField[2][blockLattice.getCellId(iX, iY+1, iZ-1)]
166 + rhoField[2][blockLattice.getCellId(iX-1, iY-1, iZ)]
167 + 2. * rhoField[2][blockLattice.getCellId(iX, iY-1, iZ)]
168 + rhoField[2][blockLattice.getCellId(iX+1, iY-1, iZ)]
169 + 2. * rhoField[2][blockLattice.getCellId(iX-1, iY, iZ)]
170 -24. * rhoField[2][blockLattice.getCellId(iX, iY, iZ)]
171 + 2. * rhoField[2][blockLattice.getCellId(iX+1, iY, iZ)]
172 + rhoField[2][blockLattice.getCellId(iX-1, iY+1, iZ)]
173 + 2. * rhoField[2][blockLattice.getCellId(iX, iY+1, iZ)]
174 + rhoField[2][blockLattice.getCellId(iX+1, iY+1, iZ)]
175 + rhoField[2][blockLattice.getCellId(iX, iY-1, iZ+1)]
176 + rhoField[2][blockLattice.getCellId(iX-1, iY, iZ+1)]
177 + 2. * rhoField[2][blockLattice.getCellId(iX, iY, iZ+1)]
178 + rhoField[2][blockLattice.getCellId(iX+1, iY, iZ+1)]
179 + rhoField[2][blockLattice.getCellId(iX, iY+1, iZ+1)]
180 );
181 }
182
183 // setting chemical potential to the respective lattices
184 blockLattice.get(iX,iY,iZ).template setField<descriptors::CHEM_POTENTIAL>(
185 term1 + term2
186 + 0.25*alpha*alpha*( (kappa2 - kappa1) * laplaceRho2
187 +(kappa2 + kappa1) * (laplaceRho3 - laplaceRho1) )
188 );
189 partnerLattice1->get(iX, iY, iZ).template setField<descriptors::CHEM_POTENTIAL>(
190 term1 - term2
191 + 0.25*alpha*alpha*( (kappa2 - kappa1) * (laplaceRho1 - laplaceRho3)
192 -(kappa2 + kappa1) * laplaceRho2 )
193 );
194 if (partners.size() > 1) {
195 partnerLattice2->get(iX, iY, iZ).template setField<descriptors::CHEM_POTENTIAL>(
196 - term1 - term2 + term3
197 + 0.25*alpha*alpha*( (kappa2 + kappa1) * laplaceRho1
198 -(kappa2 - kappa1) * laplaceRho2
199 -(kappa2 + kappa1 + 4.*kappa3) * laplaceRho3 )
200 );
201 }
202 }
203 }
204 }
205 }
206}
207
208template<typename T, typename DESCRIPTOR>
209void FreeEnergyChemicalPotentialCoupling3D<T,DESCRIPTOR>::process (
210 BlockLattice<T,DESCRIPTOR>& blockLattice)
211{
212 processSubDomain(blockLattice, x0, x1, y0, y1, z0, z1);
213}
214
216
217template<typename T, typename DESCRIPTOR>
218FreeEnergyForceCoupling3D <T,DESCRIPTOR>::FreeEnergyForceCoupling3D (
219 int x0_, int x1_, int y0_, int y1_, int z0_, int z1_,
220 std::vector<BlockStructureD<3>*> partners_)
221 : x0(x0_), x1(x1_), y0(y0_), y1(y1_), z0(z0_), z1(z1_), partners(partners_)
222{
223 this->getName() = "FreeEnergyForceCoupling3D";
224}
225
226template<typename T, typename DESCRIPTOR>
227FreeEnergyForceCoupling3D <T,DESCRIPTOR>::FreeEnergyForceCoupling3D (
228 std::vector<BlockStructureD<3>*> partners_)
229 : x0(0), x1(0), y0(0), y1(0), z0(0), z1(0), partners(partners_)
230{
231 this->getName() = "FreeEnergyForceCoupling3D";
232}
233
234template<typename T, typename DESCRIPTOR>
235void FreeEnergyForceCoupling3D<T,DESCRIPTOR>::processSubDomain (
236 BlockLattice<T,DESCRIPTOR>& blockLattice,
237 int x0_, int x1_, int y0_, int y1_, int z0_, int z1_ )
238{
239 // If partners.size() == 1: two fluid components
240 // If partners.size() == 2: three fluid components
241 BlockLattice<T,DESCRIPTOR> *partnerLattice1 = static_cast<BlockLattice<T,DESCRIPTOR> *>(partners[0]);
242 BlockLattice<T,DESCRIPTOR> *partnerLattice2 = 0;
243 if (partners.size() > 1) {
244 partnerLattice2 = static_cast<BlockLattice<T,DESCRIPTOR> *>(partners[1]);
245 }
246
247 int newX0, newX1, newY0, newY1, newZ0, newZ1;
248 if ( util::intersect ( x0, x1, y0, y1, z0, z1,
249 x0_, x1_, y0_, y1_, z0_, z1_,
250 newX0, newX1, newY0, newY1, newZ0, newZ1 ) ) {
251
252 for (int iX=newX0; iX<=newX1; ++iX) {
253 for (int iY=newY0; iY<=newY1; ++iY) {
254 for (int iZ=newZ0; iZ<=newZ1; ++iZ) {
255 T phi = blockLattice.get(iX,iY,iZ).computeRho();
256 T rho = partnerLattice1->get(iX,iY,iZ).computeRho();
257
258 T gradMuPhiX = 1./12. * ( -blockLattice.get(iX-1,iY,iZ-1).template getField<descriptors::CHEM_POTENTIAL>()
259 - blockLattice.get(iX-1,iY,iZ+1).template getField<descriptors::CHEM_POTENTIAL>()
260 - blockLattice.get(iX-1,iY-1,iZ ).template getField<descriptors::CHEM_POTENTIAL>()
261 - 2.* blockLattice.get(iX-1,iY,iZ ).template getField<descriptors::CHEM_POTENTIAL>()
262 - blockLattice.get(iX-1,iY+1,iZ ).template getField<descriptors::CHEM_POTENTIAL>()
263 + blockLattice.get(iX+1,iY-1,iZ ).template getField<descriptors::CHEM_POTENTIAL>()
264 + 2.* blockLattice.get(iX+1,iY,iZ ).template getField<descriptors::CHEM_POTENTIAL>()
265 + blockLattice.get(iX+1,iY+1,iZ ).template getField<descriptors::CHEM_POTENTIAL>()
266 + blockLattice.get(iX+1,iY,iZ-1).template getField<descriptors::CHEM_POTENTIAL>()
267 + blockLattice.get(iX+1,iY,iZ+1).template getField<descriptors::CHEM_POTENTIAL>() );
268 T gradMuPhiY = 1./12. * ( -blockLattice.get(iX-1,iY-1,iZ ).template getField<descriptors::CHEM_POTENTIAL>()
269 - blockLattice.get(iX+1,iY-1,iZ ).template getField<descriptors::CHEM_POTENTIAL>()
270 - blockLattice.get(iX,iY-1,iZ-1).template getField<descriptors::CHEM_POTENTIAL>()
271 - 2.* blockLattice.get(iX,iY-1,iZ ).template getField<descriptors::CHEM_POTENTIAL>()
272 - blockLattice.get(iX,iY-1,iZ+1).template getField<descriptors::CHEM_POTENTIAL>()
273 + blockLattice.get(iX,iY+1,iZ-1).template getField<descriptors::CHEM_POTENTIAL>()
274 + 2.* blockLattice.get(iX,iY+1,iZ ).template getField<descriptors::CHEM_POTENTIAL>()
275 + blockLattice.get(iX,iY+1,iZ+1).template getField<descriptors::CHEM_POTENTIAL>()
276 + blockLattice.get(iX-1,iY+1,iZ ).template getField<descriptors::CHEM_POTENTIAL>()
277 + blockLattice.get(iX+1,iY+1,iZ ).template getField<descriptors::CHEM_POTENTIAL>() );
278 T gradMuPhiZ = 1./12. * ( -blockLattice.get(iX,iY-1,iZ-1).template getField<descriptors::CHEM_POTENTIAL>()
279 - blockLattice.get(iX,iY+1,iZ-1).template getField<descriptors::CHEM_POTENTIAL>()
280 - blockLattice.get(iX-1,iY,iZ-1).template getField<descriptors::CHEM_POTENTIAL>()
281 - 2.* blockLattice.get(iX,iY,iZ-1).template getField<descriptors::CHEM_POTENTIAL>()
282 - blockLattice.get(iX+1,iY,iZ-1).template getField<descriptors::CHEM_POTENTIAL>()
283 + blockLattice.get(iX-1,iY,iZ+1).template getField<descriptors::CHEM_POTENTIAL>()
284 + 2.* blockLattice.get(iX,iY,iZ+1).template getField<descriptors::CHEM_POTENTIAL>()
285 + blockLattice.get(iX+1,iY,iZ+1).template getField<descriptors::CHEM_POTENTIAL>()
286 + blockLattice.get(iX,iY-1,iZ+1).template getField<descriptors::CHEM_POTENTIAL>()
287 + blockLattice.get(iX,iY+1,iZ+1).template getField<descriptors::CHEM_POTENTIAL>() );
288
289 T gradMuRhoX = 1./12. * ( -partnerLattice1->get(iX-1,iY,iZ-1).template getField<descriptors::CHEM_POTENTIAL>()
290 - partnerLattice1->get(iX-1,iY,iZ+1).template getField<descriptors::CHEM_POTENTIAL>()
291 - partnerLattice1->get(iX-1,iY-1,iZ ).template getField<descriptors::CHEM_POTENTIAL>()
292 - 2.* partnerLattice1->get(iX-1,iY,iZ ).template getField<descriptors::CHEM_POTENTIAL>()
293 - partnerLattice1->get(iX-1,iY+1,iZ ).template getField<descriptors::CHEM_POTENTIAL>()
294 + partnerLattice1->get(iX+1,iY-1,iZ ).template getField<descriptors::CHEM_POTENTIAL>()
295 + 2.* partnerLattice1->get(iX+1,iY,iZ ).template getField<descriptors::CHEM_POTENTIAL>()
296 + partnerLattice1->get(iX+1,iY+1,iZ ).template getField<descriptors::CHEM_POTENTIAL>()
297 + partnerLattice1->get(iX+1,iY,iZ-1).template getField<descriptors::CHEM_POTENTIAL>()
298 + partnerLattice1->get(iX+1,iY,iZ+1).template getField<descriptors::CHEM_POTENTIAL>() );
299 T gradMuRhoY = 1./12. * ( -partnerLattice1->get(iX-1,iY-1,iZ ).template getField<descriptors::CHEM_POTENTIAL>()
300 - partnerLattice1->get(iX+1,iY-1,iZ ).template getField<descriptors::CHEM_POTENTIAL>()
301 - partnerLattice1->get(iX,iY-1,iZ-1).template getField<descriptors::CHEM_POTENTIAL>()
302 - 2.* partnerLattice1->get(iX,iY-1,iZ ).template getField<descriptors::CHEM_POTENTIAL>()
303 - partnerLattice1->get(iX,iY-1,iZ+1).template getField<descriptors::CHEM_POTENTIAL>()
304 + partnerLattice1->get(iX,iY+1,iZ-1).template getField<descriptors::CHEM_POTENTIAL>()
305 + 2.* partnerLattice1->get(iX,iY+1,iZ ).template getField<descriptors::CHEM_POTENTIAL>()
306 + partnerLattice1->get(iX,iY+1,iZ+1).template getField<descriptors::CHEM_POTENTIAL>()
307 + partnerLattice1->get(iX-1,iY+1,iZ ).template getField<descriptors::CHEM_POTENTIAL>()
308 + partnerLattice1->get(iX+1,iY+1,iZ ).template getField<descriptors::CHEM_POTENTIAL>() );
309 T gradMuRhoZ = 1./12. * ( -partnerLattice1->get(iX,iY-1,iZ-1).template getField<descriptors::CHEM_POTENTIAL>()
310 - partnerLattice1->get(iX,iY+1,iZ-1).template getField<descriptors::CHEM_POTENTIAL>()
311 - partnerLattice1->get(iX-1,iY,iZ-1).template getField<descriptors::CHEM_POTENTIAL>()
312 - 2.* partnerLattice1->get(iX,iY,iZ-1).template getField<descriptors::CHEM_POTENTIAL>()
313 - partnerLattice1->get(iX+1,iY,iZ-1).template getField<descriptors::CHEM_POTENTIAL>()
314 + partnerLattice1->get(iX-1,iY,iZ+1).template getField<descriptors::CHEM_POTENTIAL>()
315 + 2.* partnerLattice1->get(iX,iY,iZ+1).template getField<descriptors::CHEM_POTENTIAL>()
316 + partnerLattice1->get(iX+1,iY,iZ+1).template getField<descriptors::CHEM_POTENTIAL>()
317 + partnerLattice1->get(iX,iY-1,iZ+1).template getField<descriptors::CHEM_POTENTIAL>()
318 + partnerLattice1->get(iX,iY+1,iZ+1).template getField<descriptors::CHEM_POTENTIAL>() );
319
320 T psi = 0.;
321 T gradMuPsiX = 0.;
322 T gradMuPsiY = 0.;
323 T gradMuPsiZ = 0.;
324 if (partners.size() > 1) {
325 psi = partnerLattice2->get(iX,iY,iZ).computeRho();
326 gradMuPsiX = 1./12. * ( -partnerLattice2->get(iX-1,iY,iZ-1).template getField<descriptors::CHEM_POTENTIAL>()
327 - partnerLattice2->get(iX-1,iY,iZ+1).template getField<descriptors::CHEM_POTENTIAL>()
328 - partnerLattice2->get(iX-1,iY-1,iZ ).template getField<descriptors::CHEM_POTENTIAL>()
329 - 2.* partnerLattice2->get(iX-1,iY,iZ ).template getField<descriptors::CHEM_POTENTIAL>()
330 - partnerLattice2->get(iX-1,iY+1,iZ ).template getField<descriptors::CHEM_POTENTIAL>()
331 + partnerLattice2->get(iX+1,iY-1,iZ ).template getField<descriptors::CHEM_POTENTIAL>()
332 + 2.* partnerLattice2->get(iX+1,iY,iZ ).template getField<descriptors::CHEM_POTENTIAL>()
333 + partnerLattice2->get(iX+1,iY+1,iZ ).template getField<descriptors::CHEM_POTENTIAL>()
334 + partnerLattice2->get(iX+1,iY,iZ-1).template getField<descriptors::CHEM_POTENTIAL>()
335 + partnerLattice2->get(iX+1,iY,iZ+1).template getField<descriptors::CHEM_POTENTIAL>() );
336 gradMuPsiY = 1./12. * ( -partnerLattice2->get(iX-1,iY-1,iZ ).template getField<descriptors::CHEM_POTENTIAL>()
337 - partnerLattice2->get(iX+1,iY-1,iZ ).template getField<descriptors::CHEM_POTENTIAL>()
338 - partnerLattice2->get(iX,iY-1,iZ-1).template getField<descriptors::CHEM_POTENTIAL>()
339 - 2.* partnerLattice2->get(iX,iY-1,iZ ).template getField<descriptors::CHEM_POTENTIAL>()
340 - partnerLattice2->get(iX,iY-1,iZ+1).template getField<descriptors::CHEM_POTENTIAL>()
341 + partnerLattice2->get(iX,iY+1,iZ-1).template getField<descriptors::CHEM_POTENTIAL>()
342 + 2.* partnerLattice2->get(iX,iY+1,iZ ).template getField<descriptors::CHEM_POTENTIAL>()
343 + partnerLattice2->get(iX,iY+1,iZ+1).template getField<descriptors::CHEM_POTENTIAL>()
344 + partnerLattice2->get(iX-1,iY+1,iZ ).template getField<descriptors::CHEM_POTENTIAL>()
345 + partnerLattice2->get(iX+1,iY+1,iZ ).template getField<descriptors::CHEM_POTENTIAL>() );
346 gradMuPsiZ = 1./12. * ( -partnerLattice2->get(iX,iY-1,iZ-1).template getField<descriptors::CHEM_POTENTIAL>()
347 - partnerLattice2->get(iX,iY+1,iZ-1).template getField<descriptors::CHEM_POTENTIAL>()
348 - partnerLattice2->get(iX-1,iY,iZ-1).template getField<descriptors::CHEM_POTENTIAL>()
349 - 2.* partnerLattice2->get(iX,iY,iZ-1).template getField<descriptors::CHEM_POTENTIAL>()
350 - partnerLattice2->get(iX+1,iY,iZ-1).template getField<descriptors::CHEM_POTENTIAL>()
351 + partnerLattice2->get(iX-1,iY,iZ+1).template getField<descriptors::CHEM_POTENTIAL>()
352 + 2.* partnerLattice2->get(iX,iY,iZ+1).template getField<descriptors::CHEM_POTENTIAL>()
353 + partnerLattice2->get(iX+1,iY,iZ+1).template getField<descriptors::CHEM_POTENTIAL>()
354 + partnerLattice2->get(iX,iY-1,iZ+1).template getField<descriptors::CHEM_POTENTIAL>()
355 + partnerLattice2->get(iX,iY+1,iZ+1).template getField<descriptors::CHEM_POTENTIAL>() );
356 }
357
358 T forceX = -rho*gradMuRhoX - phi*gradMuPhiX - psi*gradMuPsiX;
359 T forceY = -rho*gradMuRhoY - phi*gradMuPhiY - psi*gradMuPsiY;
360 T forceZ = -rho*gradMuRhoZ - phi*gradMuPhiZ - psi*gradMuPsiZ;
361 partnerLattice1->get(iX,iY,iZ).template setField<descriptors::FORCE>({forceX, forceY, forceZ});
362 T u[3];
363 partnerLattice1->get(iX,iY,iZ).computeU(u);
364 blockLattice.get(iX,iY,iZ).template setField<descriptors::FORCE>(u);
365 if (partners.size() > 1) {
366 partnerLattice2->get(iX,iY,iZ).template setField<descriptors::FORCE>(u);
367 }
368 }
369 }
370 }
371 }
372}
373
374template<typename T, typename DESCRIPTOR>
375void FreeEnergyForceCoupling3D<T,DESCRIPTOR>::process (
376 BlockLattice<T,DESCRIPTOR>& blockLattice)
377{
378 processSubDomain(blockLattice, x0, x1, y0, y1, z0, z1);
379}
380
381
383
384template<typename T, typename DESCRIPTOR>
385FreeEnergyInletOutletCoupling3D <T,DESCRIPTOR>::FreeEnergyInletOutletCoupling3D (
386 int x0_, int x1_, int y0_, int y1_, int z0_, int z1_,
387 std::vector<BlockStructureD<3>*> partners_)
388 : x0(x0_), x1(x1_), y0(y0_), y1(y1_), z0(z0_), z1(z1_), partners(partners_)
389{
390 this->getName() = "FreeEnergyInletOutletCoupling3D";
391}
392
393template<typename T, typename DESCRIPTOR>
394FreeEnergyInletOutletCoupling3D <T,DESCRIPTOR>::FreeEnergyInletOutletCoupling3D (
395 std::vector<BlockStructureD<3>*> partners_)
396 : x0(0), x1(0), y0(0), y1(0), z0(0), z1(0), partners(partners_)
397{
398 this->getName() = "FreeEnergyInletOutletCoupling3D";
399}
400
401template<typename T, typename DESCRIPTOR>
402void FreeEnergyInletOutletCoupling3D<T,DESCRIPTOR>::processSubDomain (
403 BlockLattice<T,DESCRIPTOR>& blockLattice,
404 int x0_, int x1_, int y0_, int y1_, int z0_, int z1_ )
405{
406 // If partners.size() == 1: two fluid components
407 // If partners.size() == 2: three fluid components
408 BlockLattice<T,DESCRIPTOR> *partnerLattice1 = static_cast<BlockLattice<T,DESCRIPTOR> *>(partners[0]);
409 BlockLattice<T,DESCRIPTOR> *partnerLattice2 = 0;
410 if (partners.size() > 1) {
411 partnerLattice2 = static_cast<BlockLattice<T,DESCRIPTOR> *>(partners[1]);
412 }
413
414 int newX0, newX1, newY0, newY1, newZ0, newZ1;
415 if ( util::intersect ( x0, x1, y0, y1, z0, z1,
416 x0_, x1_, y0_, y1_, z0_, z1_,
417 newX0, newX1, newY0, newY1, newZ0, newZ1 ) ) {
418 for (int iX=newX0; iX<=newX1; ++iX) {
419 for (int iY=newY0; iY<=newY1; ++iY) {
420 for (int iZ=newZ0; iZ<=newZ1; ++iZ) {
421 T u[DESCRIPTOR::d];
422 partnerLattice1->get(iX,iY,iZ).computeU(u);
423 blockLattice.get(iX,iY,iZ).defineU(u);
424 if (partners.size() > 1) {
425 partnerLattice2->get(iX,iY,iZ).defineU(u);
426 }
427 }
428 }
429 }
430 }
431}
432
433template<typename T, typename DESCRIPTOR>
434void FreeEnergyInletOutletCoupling3D<T,DESCRIPTOR>::process(
435 BlockLattice<T,DESCRIPTOR>& blockLattice)
436{
437 processSubDomain(blockLattice, x0, x1, y0, y1, z0, z1);
438}
439
440
442
443template<typename T, typename DESCRIPTOR>
444FreeEnergyDensityOutletCoupling3D <T,DESCRIPTOR>::FreeEnergyDensityOutletCoupling3D (
445 int x0_, int x1_, int y0_, int y1_, int z0_, int z1_, T rho_,
446 std::vector<BlockStructureD<3>*> partners_)
447 : x0(x0_), x1(x1_), y0(y0_), y1(y1_), z0(z0_), z1(z1_),
448 rho(rho_), partners(partners_)
449{
450 this->getName() = "FreeEnergyDensityOutletCoupling3D";
451}
452
453template<typename T, typename DESCRIPTOR>
454FreeEnergyDensityOutletCoupling3D <T,DESCRIPTOR>::FreeEnergyDensityOutletCoupling3D (
455 T rho_, std::vector<BlockStructureD<3>*> partners_)
456 : x0(0), x1(0), y0(0), y1(0), z0(0), z1(0), rho(rho_), partners(partners_)
457{
458 this->getName() = "FreeEnergyDensityOutletCoupling3D";
459}
460
461template<typename T, typename DESCRIPTOR>
462void FreeEnergyDensityOutletCoupling3D<T,DESCRIPTOR>::processSubDomain (
463 BlockLattice<T,DESCRIPTOR>& blockLattice,
464 int x0_, int x1_, int y0_, int y1_, int z0_, int z1_ )
465{
466 // If partners.size() == 1: two fluid components
467 // If partners.size() == 2: three fluid components
468 BlockLattice<T,DESCRIPTOR> *partnerLattice1 = static_cast<BlockLattice<T,DESCRIPTOR> *>(partners[0]);
469 BlockLattice<T,DESCRIPTOR> *partnerLattice2 = 0;
470 if (partners.size() > 1) {
471 partnerLattice2 = static_cast<BlockLattice<T,DESCRIPTOR> *>(partners[1]);
472 }
473
474 int newX0, newX1, newY0, newY1, newZ0, newZ1;
475 if ( util::intersect ( x0, x1, y0, y1, z0, z1,
476 x0_, x1_, y0_, y1_, z0_, z1_,
477 newX0, newX1, newY0, newY1, newZ0, newZ1 ) ) {
478 for (int iX=newX0; iX<=newX1; ++iX) {
479 for (int iY=newY0; iY<=newY1; ++iY) {
480 for (int iZ=newZ0; iZ<=newZ1; ++iZ) {
481
482 T rho0, phi, psi;
483 rho0 = blockLattice.get(iX,iY,iZ).computeRho();
484 phi = partnerLattice1->get(iX,iY,iZ).computeRho();
485 blockLattice.get(iX,iY,iZ).defineRho(rho);
486 partnerLattice1->get(iX,iY,iZ).defineRho(phi * rho / rho0);
487 if (partners.size() > 1) {
488 psi = partnerLattice2->get(iX,iY,iZ).computeRho();
489 partnerLattice2->get(iX,iY,iZ).defineRho(psi * rho / rho0);
490 }
491
492 }
493 }
494 }
495 }
496}
497
498template<typename T, typename DESCRIPTOR>
499void FreeEnergyDensityOutletCoupling3D<T,DESCRIPTOR>::process(
500 BlockLattice<T,DESCRIPTOR>& blockLattice)
501{
502 processSubDomain(blockLattice, x0, x1, y0, y1, z0, z1);
503}
504
505
507
508template<typename T, typename DESCRIPTOR>
509FreeEnergyChemicalPotentialGenerator3D<T,DESCRIPTOR>::FreeEnergyChemicalPotentialGenerator3D (
510 int x0_, int x1_, int y0_, int y1_, int z0_, int z1_,
511 T alpha_, T kappa1_, T kappa2_ )
512 : LatticeCouplingGenerator3D<T,DESCRIPTOR>(x0_, x1_, y0_, y1_, z0_, z1_), alpha(alpha_),
513 kappa1(kappa1_), kappa2(kappa2_), kappa3(0)
514{ }
515
516template<typename T, typename DESCRIPTOR>
517FreeEnergyChemicalPotentialGenerator3D<T,DESCRIPTOR>::FreeEnergyChemicalPotentialGenerator3D (
518 T alpha_, T kappa1_, T kappa2_ )
519 : LatticeCouplingGenerator3D<T,DESCRIPTOR>(0, 0, 0, 0, 0, 0), alpha(alpha_),
520 kappa1(kappa1_), kappa2(kappa2_), kappa3(0)
521{ }
522
523template<typename T, typename DESCRIPTOR>
524FreeEnergyChemicalPotentialGenerator3D<T,DESCRIPTOR>::FreeEnergyChemicalPotentialGenerator3D (
525 int x0_, int x1_, int y0_, int y1_, int z0_, int z1_,
526 T alpha_, T kappa1_, T kappa2_, T kappa3_ )
527 : LatticeCouplingGenerator3D<T,DESCRIPTOR>(x0_, x1_, y0_, y1_, z0_, z1_), alpha(alpha_),
528 kappa1(kappa1_), kappa2(kappa2_), kappa3(kappa3_)
529{ }
530
531template<typename T, typename DESCRIPTOR>
532FreeEnergyChemicalPotentialGenerator3D<T,DESCRIPTOR>::FreeEnergyChemicalPotentialGenerator3D (
533 T alpha_, T kappa1_, T kappa2_, T kappa3_ )
534 : LatticeCouplingGenerator3D<T,DESCRIPTOR>(0, 0, 0, 0, 0, 0), alpha(alpha_),
535 kappa1(kappa1_), kappa2(kappa2_), kappa3(kappa3_)
536{ }
537
538template<typename T, typename DESCRIPTOR>
539PostProcessor3D<T,DESCRIPTOR>* FreeEnergyChemicalPotentialGenerator3D<T,DESCRIPTOR>::generate (
540 std::vector<BlockStructureD<3>*> partners) const
541{
542 return new FreeEnergyChemicalPotentialCoupling3D<T,DESCRIPTOR>(
543 this->x0, this->x1, this->y0, this->y1, this->z0, this->z1,
544 alpha, kappa1, kappa2, kappa3, partners);
545}
546
547template<typename T, typename DESCRIPTOR>
548LatticeCouplingGenerator3D<T,DESCRIPTOR>*
553
555
556template<typename T, typename DESCRIPTOR>
557FreeEnergyForceGenerator3D<T,DESCRIPTOR>::FreeEnergyForceGenerator3D (
558 int x0_, int x1_, int y0_, int y1_, int z0_, int z1_)
559 : LatticeCouplingGenerator3D<T,DESCRIPTOR>(x0_, x1_, y0_, y1_, z0_, z1_)
560{ }
561
562template<typename T, typename DESCRIPTOR>
566
567template<typename T, typename DESCRIPTOR>
568PostProcessor3D<T,DESCRIPTOR>* FreeEnergyForceGenerator3D<T,DESCRIPTOR>::generate (
569 std::vector<BlockStructureD<3>*> partners) const
570{
571 return new FreeEnergyForceCoupling3D<T,DESCRIPTOR>(
572 this->x0, this->x1, this->y0, this->y1, this->z0, this->z1, partners);
573}
574
575template<typename T, typename DESCRIPTOR>
580
582
583template<typename T, typename DESCRIPTOR>
584FreeEnergyInletOutletGenerator3D<T,DESCRIPTOR>::FreeEnergyInletOutletGenerator3D (
585 int x0_, int x1_, int y0_, int y1_, int z0_, int z1_)
586 : LatticeCouplingGenerator3D<T,DESCRIPTOR>(x0_, x1_, y0_, y1_, z0_, z1_)
587{ }
588
589template<typename T, typename DESCRIPTOR>
593
594template<typename T, typename DESCRIPTOR>
595PostProcessor3D<T,DESCRIPTOR>* FreeEnergyInletOutletGenerator3D<T,DESCRIPTOR>::generate (
596 std::vector<BlockStructureD<3>*> partners) const
597{
598 return new FreeEnergyInletOutletCoupling3D<T,DESCRIPTOR>(
599 this->x0, this->x1, this->y0, this->y1, this->z0, this->z1, partners);
600}
601
602template<typename T, typename DESCRIPTOR>
607
609
610template<typename T, typename DESCRIPTOR>
611FreeEnergyDensityOutletGenerator3D<T,DESCRIPTOR>::FreeEnergyDensityOutletGenerator3D (
612 int x0_, int x1_, int y0_, int y1_, int z0_, int z1_, T rho_)
613 : LatticeCouplingGenerator3D<T,DESCRIPTOR>(x0_, x1_, y0_, y1_), rho(rho_)
614{ }
615
616template<typename T, typename DESCRIPTOR>
621
622template<typename T, typename DESCRIPTOR>
623PostProcessor3D<T,DESCRIPTOR>* FreeEnergyDensityOutletGenerator3D<T,DESCRIPTOR>::generate (
624 std::vector<BlockStructureD<3>*> partners) const
625{
626 return new FreeEnergyDensityOutletCoupling3D<T,DESCRIPTOR>(
627 this->x0,this->x1,this->y0,this->y1,this->z0,this->z1, rho, partners);
628}
629
630template<typename T, typename DESCRIPTOR>
635
636
637} // namespace olb
638
639#endif
olb::BlockLattice
Platform-abstracted block lattice for external access and inter-block interaction.
Definition blockLattice.h:87
olb::BlockLattice::get
Cell< T, DESCRIPTOR > get(CellID iCell)
Get Cell interface for index iCell.
Definition blockLattice.h:171
olb::BlockStructureD
Base of a regular block.
Definition blockStructure.h:54
olb::BlockStructureD::getCellId
CellID getCellId(LatticeR< D > latticeR) const
Get 1D cell ID.
Definition blockStructure.h:125
olb::FreeEnergyChemicalPotentialCoupling3D
This class calculates the chemical potential and stores it in the external field of the respective la...
Definition freeEnergyPostProcessor3D.h:41
olb::FreeEnergyChemicalPotentialCoupling3D::FreeEnergyChemicalPotentialCoupling3D
FreeEnergyChemicalPotentialCoupling3D(int x0_, int x1_, int y0_, int y1_, int z0_, int z1_, T alpha_, T kappa1_, T kappa2_, T kappa3_, std::vector< BlockStructureD< 3 > * > partners_)
Definition freeEnergyPostProcessor3D.hh:34
olb::FreeEnergyChemicalPotentialCoupling3D::process
void process(BlockLattice< T, DESCRIPTOR > &blockLattice) override
Execute post-processing step.
Definition freeEnergyPostProcessor3D.hh:209
olb::FreeEnergyChemicalPotentialCoupling3D::processSubDomain
void processSubDomain(BlockLattice< T, DESCRIPTOR > &blockLattice, int x0_, int x1_, int y0_, int y1_, int z0_, int z1_) override
Execute post-processing step on a sublattice.
Definition freeEnergyPostProcessor3D.hh:54
olb::FreeEnergyChemicalPotentialGenerator3D
Generator class for the PostProcessors calculating the chemical potential.
Definition freeEnergyPostProcessor3D.h:162
olb::FreeEnergyChemicalPotentialGenerator3D::generate
PostProcessor3D< T, DESCRIPTOR > * generate(std::vector< BlockStructureD< 3 > * > partners) const override
Definition freeEnergyPostProcessor3D.hh:539
olb::FreeEnergyChemicalPotentialGenerator3D::clone
LatticeCouplingGenerator3D< T, DESCRIPTOR > * clone() const override
Definition freeEnergyPostProcessor3D.hh:549
olb::FreeEnergyChemicalPotentialGenerator3D::FreeEnergyChemicalPotentialGenerator3D
FreeEnergyChemicalPotentialGenerator3D(int x0_, int x1_, int y0_, int y1_, int z0_, int z1_, T alpha_, T kappa1_, T kappa2_)
Two component free energy model.
Definition freeEnergyPostProcessor3D.hh:509
olb::FreeEnergyDensityOutletCoupling3D
PostProcessor for setting a constant density outlet.
Definition freeEnergyPostProcessor3D.h:133
olb::FreeEnergyDensityOutletCoupling3D::process
void process(BlockLattice< T, DESCRIPTOR > &blockLattice) override
Execute post-processing step.
Definition freeEnergyPostProcessor3D.hh:499
olb::FreeEnergyDensityOutletCoupling3D::FreeEnergyDensityOutletCoupling3D
FreeEnergyDensityOutletCoupling3D(int x0_, int x1_, int y0_, int y1_, int z0_, int z1_, T rho_, std::vector< BlockStructureD< 3 > * > partners_)
Definition freeEnergyPostProcessor3D.hh:444
olb::FreeEnergyDensityOutletCoupling3D::processSubDomain
void processSubDomain(BlockLattice< T, DESCRIPTOR > &blockLattice, int x0_, int x1_, int y0_, int y1_, int z0_, int z1_) override
Execute post-processing step on a sublattice.
Definition freeEnergyPostProcessor3D.hh:462
olb::FreeEnergyDensityOutletGenerator3D
Generator class for the PostProcessors assigning the density boundary condition at the outlet.
Definition freeEnergyPostProcessor3D.h:217
olb::FreeEnergyDensityOutletGenerator3D::clone
LatticeCouplingGenerator3D< T, DESCRIPTOR > * clone() const override
Definition freeEnergyPostProcessor3D.hh:631
olb::FreeEnergyDensityOutletGenerator3D::generate
PostProcessor3D< T, DESCRIPTOR > * generate(std::vector< BlockStructureD< 3 > * > partners) const override
Definition freeEnergyPostProcessor3D.hh:623
olb::FreeEnergyDensityOutletGenerator3D::FreeEnergyDensityOutletGenerator3D
FreeEnergyDensityOutletGenerator3D(int x0_, int x1_, int y0_, int y1_, int z0_, int z1_, T rho_)
Definition freeEnergyPostProcessor3D.hh:611
olb::FreeEnergyForceCoupling3D
PostProcessor calculating the interfacial force in the free energy model.
Definition freeEnergyPostProcessor3D.h:81
olb::FreeEnergyForceCoupling3D::FreeEnergyForceCoupling3D
FreeEnergyForceCoupling3D(int x0_, int x1_, int y0_, int y1_, int z0_, int z1_, std::vector< BlockStructureD< 3 > * > partners_)
Definition freeEnergyPostProcessor3D.hh:218
olb::FreeEnergyForceCoupling3D::process
void process(BlockLattice< T, DESCRIPTOR > &blockLattice) override
Execute post-processing step.
Definition freeEnergyPostProcessor3D.hh:375
olb::FreeEnergyForceCoupling3D::processSubDomain
void processSubDomain(BlockLattice< T, DESCRIPTOR > &blockLattice, int x0_, int x1_, int y0_, int y1_, int z0_, int z1_) override
Execute post-processing step on a sublattice.
Definition freeEnergyPostProcessor3D.hh:235
olb::FreeEnergyForceGenerator3D
Generator class for the PostProcessors calculating the interfacial force.
Definition freeEnergyPostProcessor3D.h:196
olb::FreeEnergyForceGenerator3D::generate
PostProcessor3D< T, DESCRIPTOR > * generate(std::vector< BlockStructureD< 3 > * > partners) const override
Definition freeEnergyPostProcessor3D.hh:568
olb::FreeEnergyForceGenerator3D::clone
LatticeCouplingGenerator3D< T, DESCRIPTOR > * clone() const override
Definition freeEnergyPostProcessor3D.hh:576
olb::FreeEnergyInletOutletCoupling3D
PostProcessor for assigning the velocity at inlet and outlets to lattice two and three.
Definition freeEnergyPostProcessor3D.h:106
olb::FreeEnergyInletOutletCoupling3D::processSubDomain
void processSubDomain(BlockLattice< T, DESCRIPTOR > &blockLattice, int x0_, int x1_, int y0_, int y1_, int z0_, int z1_) override
Execute post-processing step on a sublattice.
Definition freeEnergyPostProcessor3D.hh:402
olb::FreeEnergyInletOutletCoupling3D::process
void process(BlockLattice< T, DESCRIPTOR > &blockLattice) override
Execute post-processing step.
Definition freeEnergyPostProcessor3D.hh:434
olb::FreeEnergyInletOutletCoupling3D::FreeEnergyInletOutletCoupling3D
FreeEnergyInletOutletCoupling3D(int x0_, int x1_, int y0_, int y1_, int z0_, int z1_, std::vector< BlockStructureD< 3 > * > partners_)
Definition freeEnergyPostProcessor3D.hh:385
olb::FreeEnergyInletOutletGenerator3D
Generator class for the PostProcessors assigning the velocity at the outlet to lattice two and three.
Definition freeEnergyPostProcessor3D.h:206
olb::FreeEnergyInletOutletGenerator3D::clone
LatticeCouplingGenerator3D< T, DESCRIPTOR > * clone() const override
Definition freeEnergyPostProcessor3D.hh:603
olb::FreeEnergyInletOutletGenerator3D::generate
PostProcessor3D< T, DESCRIPTOR > * generate(std::vector< BlockStructureD< 3 > * > partners) const override
Definition freeEnergyPostProcessor3D.hh:595
olb::PostProcessor3D::getName
std::string & getName()
read and write access to name
Definition postProcessing.hh:192
olb::util::intersect
bool intersect(int x0, int x1, int y0, int y1, int x0_, int x1_, int y0_, int y1_, int &newX0, int &newX1, int &newY0, int &newY1)
Definition util.h:89
olb
Top level namespace for all of OpenLB.
Definition boundaryPostProcessors2D.h:34
Generated on Tue Mar 5 2024 09:57:47 for OpenLB by doxygen 1.10.0