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/freeEnergyPostProcessor2D.hh | 511 ++++++++++++++++++++++++++++++
1 file changed, 511 insertions(+)
create mode 100644 src/dynamics/freeEnergyPostProcessor2D.hh
(limited to 'src/dynamics/freeEnergyPostProcessor2D.hh')
diff --git a/src/dynamics/freeEnergyPostProcessor2D.hh b/src/dynamics/freeEnergyPostProcessor2D.hh
new file mode 100644
index 0000000..f6aa19f
--- /dev/null
+++ b/src/dynamics/freeEnergyPostProcessor2D.hh
@@ -0,0 +1,511 @@
+/* This file is part of the OpenLB library
+ *
+ * Copyright (C) 2018 Robin Trunk, Sam Avis
+ * 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.
+*/
+
+#ifndef FREE_ENERGY_POST_PROCESSOR_2D_HH
+#define FREE_ENERGY_POST_PROCESSOR_2D_HH
+
+#include "freeEnergyPostProcessor2D.h"
+#include "core/blockLattice2D.h"
+
+namespace olb {
+
+//////// FreeEnergyChemicalPotentialCoupling2D ///////////////////////////////////
+
+template
+FreeEnergyChemicalPotentialCoupling2D ::FreeEnergyChemicalPotentialCoupling2D (
+ int x0_, int x1_, int y0_, int y1_, T alpha_, T kappa1_, T kappa2_, T kappa3_,
+ std::vector partners_)
+ : x0(x0_), x1(x1_), y0(y0_), y1(y1_), alpha(alpha_), kappa1(kappa1_),
+ kappa2(kappa2_), kappa3(kappa3_), partners(partners_)
+{ }
+
+template
+FreeEnergyChemicalPotentialCoupling2D ::FreeEnergyChemicalPotentialCoupling2D (
+ T alpha_, T kappa1_, T kappa2_, T kappa3_, std::vector partners_)
+ : x0(0), x1(0), y0(0), y1(0), alpha(alpha_), kappa1(kappa1_), kappa2(kappa2_),
+ kappa3(kappa3_), partners(partners_)
+{ }
+
+template
+void FreeEnergyChemicalPotentialCoupling2D::processSubDomain (
+ BlockLattice2D& blockLattice, int x0_, int x1_, int y0_, int y1_ )
+{
+ // If partners.size() == 1: two fluid components
+ // If partners.size() == 2: three fluid components
+ BlockLattice2D *partnerLattice1 = dynamic_cast *>(partners[0]);
+ BlockLattice2D *partnerLattice2 = 0;
+ if (partners.size() > 1) {
+ partnerLattice2 = dynamic_cast *>(partners[1]);
+ }
+
+ int newX0, newX1, newY0, newY1;
+ if ( util::intersect ( x0, x1, y0, y1,
+ x0_, x1_, y0_, y1_,
+ newX0, newX1, newY0, newY1 ) ) {
+ int nx = newX1-newX0+3; // include a one-cell boundary
+ int ny = newY1-newY0+3; // include a one-cell boundary
+ int offsetX = newX0-1;
+ int offsetY = newY0-1;
+
+ // compute the density fields for each lattice
+ BlockData2D rhoField1(nx, ny);
+ BlockData2D rhoField2(nx, ny);
+ BlockData2D rhoField3(nx, ny);
+ for (int iX=newX0-1; iX<=newX1+1; ++iX)
+ for (int iY=newY0-1; iY<=newY1+1; ++iY)
+ rhoField1.get(iX-offsetX, iY-offsetY) = blockLattice.get(iX,iY).computeRho();
+ for (int iX=newX0-1; iX<=newX1+1; ++iX)
+ for (int iY=newY0-1; iY<=newY1+1; ++iY)
+ rhoField2.get(iX-offsetX, iY-offsetY) = partnerLattice1->get(iX,iY).computeRho();
+ if (partners.size() > 1) {
+ for (int iX=newX0-1; iX<=newX1+1; ++iX)
+ for (int iY=newY0-1; iY<=newY1+1; ++iY)
+ rhoField3.get(iX-offsetX, iY-offsetY) = partnerLattice2->get(iX,iY).computeRho();
+ }
+
+ // calculate chemical potential
+ for (int iX=newX0; iX<=newX1; ++iX) {
+ for (int iY=newY0; iY<=newY1; ++iY) {
+ T densitySum = rhoField1.get(iX-offsetX, iY-offsetY)
+ + rhoField2.get(iX-offsetX, iY-offsetY);
+ T densityDifference = rhoField1.get(iX-offsetX, iY-offsetY)
+ - rhoField2.get(iX-offsetX, iY-offsetY);
+ if (partners.size() > 1) {
+ densitySum -= rhoField3.get(iX-offsetX, iY-offsetY);
+ densityDifference -= rhoField3.get(iX-offsetX, iY-offsetY);
+ }
+ T term1 = 0.125 * kappa1 * (densitySum)
+ * (densitySum-1.) * (densitySum-2.);
+ T term2 = 0.125 * kappa2 * (densityDifference)
+ * (densityDifference-1.) * (densityDifference-2.);
+ T term3 = 0.;
+ if (partners.size() > 1) {
+ T rho3 = rhoField3.get(iX-offsetX, iY-offsetY);
+ term3 = kappa3 * rho3 * (rho3 - 1.) * (2.*rho3 - 1.);
+ }
+
+ T laplaceRho1 = 0.25 * (
+ rhoField1.get(iX-offsetX-1, iY-offsetY-1)
+ + 2. * rhoField1.get(iX-offsetX, iY-offsetY-1)
+ + rhoField1.get(iX-offsetX+1, iY-offsetY-1)
+ + 2. * rhoField1.get(iX-offsetX-1, iY-offsetY)
+ -12. * rhoField1.get(iX-offsetX, iY-offsetY)
+ + 2. * rhoField1.get(iX-offsetX+1, iY-offsetY)
+ + rhoField1.get(iX-offsetX-1, iY-offsetY+1)
+ + 2. * rhoField1.get(iX-offsetX, iY-offsetY+1)
+ + rhoField1.get(iX-offsetX+1, iY-offsetY+1)
+ );
+ T laplaceRho2 = 0.25 * (
+ rhoField2.get(iX-offsetX-1, iY-offsetY-1)
+ + 2. * rhoField2.get(iX-offsetX, iY-offsetY-1)
+ + rhoField2.get(iX-offsetX+1, iY-offsetY-1)
+ + 2. * rhoField2.get(iX-offsetX-1, iY-offsetY)
+ -12. * rhoField2.get(iX-offsetX, iY-offsetY)
+ + 2. * rhoField2.get(iX-offsetX+1, iY-offsetY)
+ + rhoField2.get(iX-offsetX-1, iY-offsetY+1)
+ + 2. * rhoField2.get(iX-offsetX, iY-offsetY+1)
+ + rhoField2.get(iX-offsetX+1, iY-offsetY+1)
+ );
+ T laplaceRho3 = 0.;
+ if (partners.size() > 1) {
+ laplaceRho3 = 0.25 * (
+ rhoField3.get(iX-offsetX-1, iY-offsetY-1)
+ + 2. * rhoField3.get(iX-offsetX, iY-offsetY-1)
+ + rhoField3.get(iX-offsetX+1, iY-offsetY-1)
+ + 2. * rhoField3.get(iX-offsetX-1, iY-offsetY)
+ -12. * rhoField3.get(iX-offsetX, iY-offsetY)
+ + 2. * rhoField3.get(iX-offsetX+1, iY-offsetY)
+ + rhoField3.get(iX-offsetX-1, iY-offsetY+1)
+ + 2. * rhoField3.get(iX-offsetX, iY-offsetY+1)
+ + rhoField3.get(iX-offsetX+1, iY-offsetY+1)
+ );
+ }
+
+ // setting chemical potential to the respective lattices
+ blockLattice.get(iX, iY).template setField(term1 + term2
+ + 0.25*alpha*alpha*( (kappa2 - kappa1) * laplaceRho2
+ +(kappa2 + kappa1) * (laplaceRho3 - laplaceRho1) ));
+ partnerLattice1->get(iX, iY).template setField(term1 - term2
+ + 0.25*alpha*alpha*( (kappa2 - kappa1) * (laplaceRho1 - laplaceRho3)
+ -(kappa2 + kappa1) * laplaceRho2 ));
+ if (partners.size() > 1) {
+ partnerLattice2->get(iX, iY).template setField(- term1 - term2 + term3
+ + 0.25*alpha*alpha*( (kappa2 + kappa1) * laplaceRho1
+ -(kappa2 - kappa1) * laplaceRho2
+ -(kappa2 + kappa1 + 4.*kappa3) * laplaceRho3 ));
+ }
+ }
+ }
+
+ }
+}
+
+template
+void FreeEnergyChemicalPotentialCoupling2D::process (
+ BlockLattice2D& blockLattice)
+{
+ processSubDomain(blockLattice, x0, x1, y0, y1);
+}
+
+
+//////// FreeEnergyForceCoupling2D ///////////////////////////////////
+
+template
+FreeEnergyForceCoupling2D ::FreeEnergyForceCoupling2D (
+ int x0_, int x1_, int y0_, int y1_,
+ std::vector partners_)
+ : x0(x0_), x1(x1_), y0(y0_), y1(y1_), partners(partners_)
+{ }
+
+template
+FreeEnergyForceCoupling2D ::FreeEnergyForceCoupling2D (
+ std::vector partners_)
+ : x0(0), x1(0), y0(0), y1(0), partners(partners_)
+{ }
+
+template
+void FreeEnergyForceCoupling2D::processSubDomain (
+ BlockLattice2D& blockLattice, int x0_, int x1_, int y0_, int y1_ )
+{
+ // If partners.size() == 1: two fluid components
+ // If partners.size() == 2: three fluid components
+ BlockLattice2D *partnerLattice1 = dynamic_cast *>(partners[0]);
+ BlockLattice2D *partnerLattice2 = 0;
+ if (partners.size() > 1) {
+ partnerLattice2 = dynamic_cast *>(partners[1]);
+ }
+
+ int newX0, newX1, newY0, newY1;
+ if ( util::intersect ( x0, x1, y0, y1,
+ x0_, x1_, y0_, y1_,
+ newX0, newX1, newY0, newY1 ) ) {
+
+ for (int iX=newX0; iX<=newX1; ++iX) {
+ for (int iY=newY0; iY<=newY1; ++iY) {
+ T phi = blockLattice.get(iX,iY).computeRho();
+ T rho = partnerLattice1->get(iX,iY).computeRho();
+ T gradMuPhiX = 1./12. * ( -blockLattice.get(iX-1,iY-1).template getField()
+ - 4.* blockLattice.get(iX-1,iY ).template getField()
+ - blockLattice.get(iX-1,iY+1).template getField()
+ + blockLattice.get(iX+1,iY-1).template getField()
+ + 4.* blockLattice.get(iX+1,iY ).template getField()
+ + blockLattice.get(iX+1,iY+1).template getField() );
+ T gradMuPhiY = 1./12. * ( -blockLattice.get(iX-1,iY-1).template getField()
+ - 4.* blockLattice.get(iX ,iY-1).template getField()
+ - blockLattice.get(iX+1,iY-1).template getField()
+ + blockLattice.get(iX-1,iY+1).template getField()
+ + 4.* blockLattice.get(iX ,iY+1).template getField()
+ + blockLattice.get(iX+1,iY+1).template getField() );
+ T gradMuRhoX = 1./12. * ( -partnerLattice1->get(iX-1,iY-1).template getField()
+ - 4.* partnerLattice1->get(iX-1,iY ).template getField()
+ - partnerLattice1->get(iX-1,iY+1).template getField()
+ + partnerLattice1->get(iX+1,iY-1).template getField()
+ + 4.* partnerLattice1->get(iX+1,iY ).template getField()
+ + partnerLattice1->get(iX+1,iY+1).template getField() );
+ T gradMuRhoY = 1./12. * ( -partnerLattice1->get(iX-1,iY-1).template getField()
+ - 4.* partnerLattice1->get(iX ,iY-1).template getField()
+ - partnerLattice1->get(iX+1,iY-1).template getField()
+ + partnerLattice1->get(iX-1,iY+1).template getField()
+ + 4.* partnerLattice1->get(iX ,iY+1).template getField()
+ + partnerLattice1->get(iX+1,iY+1).template getField() );
+ T psi = 0.;
+ T gradMuPsiX = 0.;
+ T gradMuPsiY = 0.;
+ if (partners.size() > 1) {
+ psi = partnerLattice2->get(iX,iY).computeRho();
+ gradMuPsiX = 1./12. * ( -partnerLattice2->get(iX-1,iY-1).template getField()
+ - 4.* partnerLattice2->get(iX-1,iY ).template getField()
+ - partnerLattice2->get(iX-1,iY+1).template getField()
+ + partnerLattice2->get(iX+1,iY-1).template getField()
+ + 4.* partnerLattice2->get(iX+1,iY ).template getField()
+ + partnerLattice2->get(iX+1,iY+1).template getField() );
+ gradMuPsiY = 1./12. * ( -partnerLattice2->get(iX-1,iY-1).template getField()
+ - 4.* partnerLattice2->get(iX ,iY-1).template getField()
+ - partnerLattice2->get(iX+1,iY-1).template getField()
+ + partnerLattice2->get(iX-1,iY+1).template getField()
+ + 4.* partnerLattice2->get(iX ,iY+1).template getField()
+ + partnerLattice2->get(iX+1,iY+1).template getField() );
+ }
+
+ T forceX = -rho*gradMuRhoX - phi*gradMuPhiX - psi*gradMuPsiX;
+ T forceY = -rho*gradMuRhoY - phi*gradMuPhiY - psi*gradMuPsiY;
+ partnerLattice1->get(iX, iY).template setField({forceX, forceY});
+ T u[2];
+ partnerLattice1->get(iX,iY).computeU(u);
+ blockLattice.get(iX, iY).template setField(u);
+ if (partners.size() > 1) {
+ partnerLattice2->get(iX, iY).template setField(u);
+ }
+ }
+ }
+ }
+}
+
+template
+void FreeEnergyForceCoupling2D::process(
+ BlockLattice2D& blockLattice)
+{
+ processSubDomain(blockLattice, x0, x1, y0, y1);
+}
+
+
+//////// FreeEnergyInletOutletCoupling2D ///////////////////////////////////
+
+template
+FreeEnergyInletOutletCoupling2D ::FreeEnergyInletOutletCoupling2D (
+ int x0_, int x1_, int y0_, int y1_, std::vector partners_)
+ : x0(x0_), x1(x1_), y0(y0_), y1(y1_), partners(partners_)
+{ }
+
+template
+FreeEnergyInletOutletCoupling2D ::FreeEnergyInletOutletCoupling2D (
+ std::vector partners_)
+ : x0(0), x1(0), y0(0), y1(0), partners(partners_)
+{ }
+
+template
+void FreeEnergyInletOutletCoupling2D::processSubDomain (
+ BlockLattice2D& blockLattice, int x0_, int x1_, int y0_, int y1_ )
+{
+ // If partners.size() == 1: two fluid components
+ // If partners.size() == 2: three fluid components
+ BlockLattice2D *partnerLattice1 = dynamic_cast *>(partners[0]);
+ BlockLattice2D *partnerLattice2 = 0;
+ if (partners.size() > 1) {
+ partnerLattice2 = dynamic_cast *>(partners[1]);
+ }
+
+ int newX0, newX1, newY0, newY1;
+ if ( util::intersect ( x0, x1, y0, y1,
+ x0_, x1_, y0_, y1_,
+ newX0, newX1, newY0, newY1 ) ) {
+ for (int iX=newX0; iX<=newX1; ++iX) {
+ for (int iY=newY0; iY<=newY1; ++iY) {
+
+ T u[DESCRIPTOR::d];
+ partnerLattice1->get(iX,iY).computeU(u);
+ blockLattice.get(iX,iY).defineU(u);
+ if (partners.size() > 1) {
+ partnerLattice2->get(iX,iY).defineU(u);
+ }
+ }
+ }
+ }
+}
+
+template
+void FreeEnergyInletOutletCoupling2D::process(
+ BlockLattice2D& blockLattice)
+{
+ processSubDomain(blockLattice, x0, x1, y0, y1);
+}
+
+
+//////// FreeEnergyDensityOutletCoupling2D ///////////////////////////////////
+
+template
+FreeEnergyDensityOutletCoupling2D ::FreeEnergyDensityOutletCoupling2D (
+ int x0_, int x1_, int y0_, int y1_, T rho_,
+ std::vector partners_)
+ : x0(x0_), x1(x1_), y0(y0_), y1(y1_), rho(rho_), partners(partners_)
+{ }
+
+template
+FreeEnergyDensityOutletCoupling2D ::FreeEnergyDensityOutletCoupling2D (
+ T rho_, std::vector partners_)
+ : x0(0), x1(0), y0(0), y1(0), rho(rho_), partners(partners_)
+{ }
+
+template
+void FreeEnergyDensityOutletCoupling2D::processSubDomain (
+ BlockLattice2D& blockLattice, int x0_, int x1_, int y0_, int y1_ )
+{
+ // If partners.size() == 1: two fluid components
+ // If partners.size() == 2: three fluid components
+ BlockLattice2D *partnerLattice1 = dynamic_cast *>(partners[0]);
+ BlockLattice2D *partnerLattice2 = 0;
+ if (partners.size() > 1) {
+ partnerLattice2 = dynamic_cast *>(partners[1]);
+ }
+
+ int newX0, newX1, newY0, newY1;
+ if ( util::intersect ( x0, x1, y0, y1,
+ x0_, x1_, y0_, y1_,
+ newX0, newX1, newY0, newY1 ) ) {
+ for (int iX=newX0; iX<=newX1; ++iX) {
+ for (int iY=newY0; iY<=newY1; ++iY) {
+
+ T rho0, phi, psi;
+ rho0 = blockLattice.get(iX,iY).computeRho();
+ phi = partnerLattice1->get(iX,iY).computeRho();
+ blockLattice.get(iX,iY).defineRho(rho);
+ partnerLattice1->get(iX,iY).defineRho(phi * rho / rho0);
+ if (partners.size() > 1) {
+ psi = partnerLattice2->get(iX,iY).computeRho();
+ partnerLattice2->get(iX,iY).defineRho(psi * rho / rho0);
+ }
+
+ }
+ }
+ }
+}
+
+template
+void FreeEnergyDensityOutletCoupling2D::process(
+ BlockLattice2D& blockLattice)
+{
+ processSubDomain(blockLattice, x0, x1, y0, y1);
+}
+
+
+//////// FreeEnergyChemicalPotentialGenerator2D ///////////////////////////////////
+
+template
+FreeEnergyChemicalPotentialGenerator2D::FreeEnergyChemicalPotentialGenerator2D (
+ int x0_, int x1_, int y0_, int y1_, T alpha_, T kappa1_, T kappa2_)
+ : LatticeCouplingGenerator2D(x0_, x1_, y0_, y1_), alpha(alpha_),
+ kappa1(kappa1_), kappa2(kappa2_), kappa3(0)
+{ }
+
+template
+FreeEnergyChemicalPotentialGenerator2D::FreeEnergyChemicalPotentialGenerator2D (
+ T alpha_, T kappa1_, T kappa2_ )
+ : LatticeCouplingGenerator2D(0, 0, 0, 0), alpha(alpha_),
+ kappa1(kappa1_), kappa2(kappa2_), kappa3(0)
+{ }
+
+template
+FreeEnergyChemicalPotentialGenerator2D::FreeEnergyChemicalPotentialGenerator2D (
+ int x0_, int x1_, int y0_, int y1_, T alpha_, T kappa1_, T kappa2_, T kappa3_)
+ : LatticeCouplingGenerator2D(x0_, x1_, y0_, y1_), alpha(alpha_),
+ kappa1(kappa1_), kappa2(kappa2_), kappa3(kappa3_)
+{ }
+
+template
+FreeEnergyChemicalPotentialGenerator2D::FreeEnergyChemicalPotentialGenerator2D (
+ T alpha_, T kappa1_, T kappa2_, T kappa3_ )
+ : LatticeCouplingGenerator2D(0, 0, 0, 0), alpha(alpha_),
+ kappa1(kappa1_), kappa2(kappa2_), kappa3(kappa3_)
+{ }
+
+template
+PostProcessor2D* FreeEnergyChemicalPotentialGenerator2D::generate (
+ std::vector partners) const
+{
+ return new FreeEnergyChemicalPotentialCoupling2D(
+ this->x0,this->x1,this->y0,this->y1, alpha, kappa1, kappa2, kappa3, partners);
+}
+
+template
+LatticeCouplingGenerator2D*
+FreeEnergyChemicalPotentialGenerator2D::clone() const
+{
+ return new FreeEnergyChemicalPotentialGenerator2D(*this);
+}
+
+//////// FreeEnergyForceGenerator2D ///////////////////////////////////
+
+template
+FreeEnergyForceGenerator2D::FreeEnergyForceGenerator2D (
+ int x0_, int x1_, int y0_, int y1_)
+ : LatticeCouplingGenerator2D(x0_, x1_, y0_, y1_)
+{ }
+
+template
+FreeEnergyForceGenerator2D::FreeEnergyForceGenerator2D ( )
+ : LatticeCouplingGenerator2D(0, 0, 0, 0)
+{ }
+
+template
+PostProcessor2D* FreeEnergyForceGenerator2D::generate (
+ std::vector partners) const
+{
+ return new FreeEnergyForceCoupling2D(
+ this->x0,this->x1,this->y0,this->y1, partners);
+}
+
+template
+LatticeCouplingGenerator2D* FreeEnergyForceGenerator2D::clone() const
+{
+ return new FreeEnergyForceGenerator2D(*this);
+}
+
+//////// FreeEnergyInletOutletGenerator2D ///////////////////////////////////
+
+template
+FreeEnergyInletOutletGenerator2D::FreeEnergyInletOutletGenerator2D (
+ int x0_, int x1_, int y0_, int y1_)
+ : LatticeCouplingGenerator2D(x0_, x1_, y0_, y1_)
+{ }
+
+template
+FreeEnergyInletOutletGenerator2D::FreeEnergyInletOutletGenerator2D ( )
+ : LatticeCouplingGenerator2D(0, 0, 0, 0)
+{ }
+
+template
+PostProcessor2D* FreeEnergyInletOutletGenerator2D::generate (
+ std::vector partners) const
+{
+ return new FreeEnergyInletOutletCoupling2D(
+ this->x0,this->x1,this->y0,this->y1, partners);
+}
+
+template
+LatticeCouplingGenerator2D* FreeEnergyInletOutletGenerator2D::clone() const
+{
+ return new FreeEnergyInletOutletGenerator2D(*this);
+}
+
+//////// FreeEnergyDensityOutletGenerator2D ///////////////////////////////////
+
+template
+FreeEnergyDensityOutletGenerator2D::FreeEnergyDensityOutletGenerator2D (
+ int x0_, int x1_, int y0_, int y1_, T rho_)
+ : LatticeCouplingGenerator2D(x0_, x1_, y0_, y1_), rho(rho_)
+{ }
+
+template
+FreeEnergyDensityOutletGenerator2D::FreeEnergyDensityOutletGenerator2D (
+ T rho_)
+ : LatticeCouplingGenerator2D(0, 0, 0, 0), rho(rho_)
+{ }
+
+template
+PostProcessor2D* FreeEnergyDensityOutletGenerator2D::generate (
+ std::vector partners) const
+{
+ return new FreeEnergyDensityOutletCoupling2D(
+ this->x0,this->x1,this->y0,this->y1, rho, partners);
+}
+
+template
+LatticeCouplingGenerator2D* FreeEnergyDensityOutletGenerator2D::clone() const
+{
+ return new FreeEnergyDensityOutletGenerator2D(*this);
+}
+
+
+} // namespace olb
+
+#endif
--
cgit v1.2.3