From 26766ef69a8e8b16d8add25d7a463a443d9f1084 Mon Sep 17 00:00:00 2001 From: James McClure Date: Sat, 25 Dec 2021 17:58:53 -0500 Subject: [PATCH] update membrane structure --- analysis/Membrane.cpp | 151 +++++++++++------------------------------- analysis/Membrane.h | 7 +- common/ScaLBL.cpp | 6 +- common/ScaLBL.h | 8 +++ 4 files changed, 52 insertions(+), 120 deletions(-) diff --git a/analysis/Membrane.cpp b/analysis/Membrane.cpp index b387dc38..a901e5e9 100644 --- a/analysis/Membrane.cpp +++ b/analysis/Membrane.cpp @@ -12,7 +12,7 @@ Membrane::~Membrane() { } -int Membrane::Create(DoubleArray &Distance, IntArray &Map, int *neighborList, int *membrane, int Np){ +int Membrane::Create(DoubleArray &Distance, IntArray &Map, int *neighborList, int Np){ int mlink = 0; int n, idx, neighbor; double dist, locdist; @@ -32,7 +32,6 @@ int Membrane::Create(DoubleArray &Distance, IntArray &Map, int *neighborList, in else if (!(idx<0)){ - int neighbor; // cycle through the neighbors of lattice site idx neighbor=Map(i-1,j,k); dist=Distance(i-1,j,k); if (dist*locdist < 0.0){ @@ -155,7 +154,8 @@ int Membrane::Create(DoubleArray &Distance, IntArray &Map, int *neighborList, in } /* allocate memory */ - membrane = new int [2*mlink]; + membraneLinks = new int [2*mlink]; + membraneDist = new double [2*mlink]; /* construct the membrane*/ mlink = 0; @@ -167,167 +167,94 @@ int Membrane::Create(DoubleArray &Distance, IntArray &Map, int *neighborList, in locdist=Distance(i,j,k); else if (!(idx<0)){ - - int neighbor; // cycle through the neighbors of lattice site idx - neighbor=Map(i-1,j,k); - dist=Distance(i-1,j,k); - if (dist*locdist < 0.0){ - neighborList[idx]=idx + 2*Np; - //membrane[2*mlink] = idx + 2*Np; - //membrane[2*mlink+1] = neighbor + 1*Np; - //mlink++; - } neighbor=Map(i+1,j,k); dist=Distance(i+1,j,k); if (dist*locdist < 0.0){ - neighborList[Np+idx] = idx + 1*Np; - membrane[2*mlink] = idx + 1*Np; - membrane[2*mlink+1] = neighbor + 2*Np; + membraneLinks[2*mlink] = idx + 1*Np; + membraneLinks[2*mlink+1] = neighbor + 2*Np; + membraneDist[2*mlink] = locdist; + membraneDist[2*mlink+1] = dist; mlink++; } - neighbor=Map(i,j-1,k); - dist=Distance(i,j-1,k); - if (dist*locdist < 0.0){ - neighborList[2*Np+idx]=idx + 4*Np; - //membrane[2*mlink] = idx + 4*Np; - //membrane[2*mlink+1] = neighbor + 3*Np; - //mlink++; - } - neighbor=Map(i,j+1,k); dist=Distance(i,j+1,k); if (dist*locdist < 0.0){ - neighborList[3*Np+idx]=idx + 3*Np; - membrane[2*mlink] = idx + 3*Np; - membrane[2*mlink+1] = neighbor + 4*Np; + membraneLinks[2*mlink] = idx + 3*Np; + membraneLinks[2*mlink+1] = neighbor + 4*Np; + membraneDist[2*mlink] = locdist; + membraneDist[2*mlink+1] = dist; mlink++; } - neighbor=Map(i,j,k-1); - dist=Distance(i,j,k-1); - if (dist*locdist < 0.0){ - neighborList[4*Np+idx]=idx + 6*Np; - //membrane[2*mlink] = idx + 6*Np; - //membrane[2*mlink+1] = neighbor + 5*Np; - //mlink++; - } - neighbor=Map(i,j,k+1); dist=Distance(i,j,k+1); if (dist*locdist < 0.0){ - neighborList[5*Np+idx]=idx + 5*Np; - membrane[2*mlink] = idx + 5*Np; - membrane[2*mlink+1] = neighbor + 6*Np; + membraneLinks[2*mlink] = idx + 5*Np; + membraneLinks[2*mlink+1] = neighbor + 6*Np; + membraneDist[2*mlink] = locdist; + membraneDist[2*mlink+1] = dist; mlink++; } - neighbor=Map(i-1,j-1,k); - dist=Distance(i-1,j-1,k); - if (dist*locdist < 0.0){ - neighborList[6*Np+idx]=idx + 8*Np; - //membrane[2*mlink] = idx + 8*Np; - //membrane[2*mlink+1] = neighbor + 7*Np - //mlink++; - } - neighbor=Map(i+1,j+1,k); dist=Distance(i+1,j+1,k); if (dist*locdist < 0.0){ - neighborList[7*Np+idx]=idx + 7*Np; - membrane[2*mlink] = idx + 7*Np; - membrane[2*mlink+1] = neighbor+8*Np; + membraneLinks[2*mlink] = idx + 7*Np; + membraneLinks[2*mlink+1] = neighbor+8*Np; + membraneDist[2*mlink] = locdist; + membraneDist[2*mlink+1] = dist; mlink++; } - neighbor=Map(i-1,j+1,k); - dist=Distance(i-1,j+1,k); - if (dist*locdist < 0.0){ - neighborList[8*Np+idx]=idx + 10*Np; - //membrane[2*mlink] = idx + 10*Np; - //membrane[2*mlink+1] = neighbor + 9*Np; - //mlink++; - } - neighbor=Map(i+1,j-1,k); dist=Distance(i+1,j-1,k); if (dist*locdist < 0.0){ - neighborList[9*Np+idx]=idx + 9*Np; - membrane[2*mlink] = idx + 9*Np; - membrane[2*mlink+1] = neighbor + 10*Np; + membraneLinks[2*mlink] = idx + 9*Np; + membraneLinks[2*mlink+1] = neighbor + 10*Np; + membraneDist[2*mlink] = locdist; + membraneDist[2*mlink+1] = dist; mlink++; } - neighbor=Map(i-1,j,k-1); - dist=Distance(i-1,j,k-1); - if (dist*locdist < 0.0){ - neighborList[10*Np+idx]=idx + 12*Np; - //membrane[2*mlink] = idx + 12*Np; - //membrane[2*mlink+1] = neighbor + 11*Np;; - //mlink++; - } - neighbor=Map(i+1,j,k+1); dist=Distance(i+1,j,k+1); if (dist*locdist < 0.0){ - neighborList[11*Np+idx]=idx + 11*Np; - membrane[2*mlink] = idx + 11*Np; - membrane[2*mlink+1] = neighbor + 12*Np; + membraneLinks[2*mlink] = idx + 11*Np; + membraneLinks[2*mlink+1] = neighbor + 12*Np; + membraneDist[2*mlink] = locdist; + membraneDist[2*mlink+1] = dist; mlink++; } - neighbor=Map(i-1,j,k+1); - dist=Distance(i-1,j,k+1); - if (dist*locdist < 0.0){ - neighborList[12*Np+idx]=idx + 14*Np; - //membrane[2*mlink] = idx + 14*Np; - //membrane[2*mlink+1] = neighbor + 13*Np; - //mlink++; - } - neighbor=Map(i+1,j,k-1); dist=Distance(i+1,j,k-1); if (dist*locdist < 0.0){ - neighborList[13*Np+idx]=idx + 13*Np; - membrane[2*mlink] = idx + 13*Np; - membrane[2*mlink+1] = neighbor + 14*Np; + membraneLinks[2*mlink] = idx + 13*Np; + membraneLinks[2*mlink+1] = neighbor + 14*Np; + membraneDist[2*mlink] = locdist; + membraneDist[2*mlink+1] = dist; mlink++; } - neighbor=Map(i,j-1,k-1); - dist=Distance(i,j-1,k-1); - if (dist*locdist < 0.0){ - neighborList[14*Np+idx]=idx + 16*Np; - //membrane[2*mlink] = idx + 16*Np; - //membrane[2*mlink+1] = neighbor + 15*Np; - //mlink++; - } - neighbor=Map(i,j+1,k+1); dist=Distance(i,j+1,k+1); if (dist*locdist < 0.0){ - neighborList[15*Np+idx]=idx + 15*Np; - membrane[2*mlink] = idx + 15*Np; - membrane[2*mlink+1] =neighbor + 16*Np; + membraneLinks[2*mlink] = idx + 15*Np; + membraneLinks[2*mlink+1] =neighbor + 16*Np; + membraneDist[2*mlink] = locdist; + membraneDist[2*mlink+1] = dist; mlink++; } - neighbor=Map(i,j-1,k+1); - dist=Distance(i,j-1,k+1); - if (dist*locdist < 0.0){ - neighborList[16*Np+idx]=idx + 18*Np; - //membrane[2*mlink] = idx + 18*Np; - //membrane[2*mlink+1] = neighbor + 17*Np; - //mlink++; - } - neighbor=Map(i,j+1,k-1); dist=Distance(i,j+1,k-1); if (dist*locdist < 0.0){ - neighborList[17*Np+idx]=idx + 17*Np; - membrane[2*mlink] = idx + 17*Np; - membrane[2*mlink+1] = neighbor + 18*Np; + membraneLinks[2*mlink] = idx + 17*Np; + membraneLinks[2*mlink+1] = neighbor + 18*Np; + membraneDist[2*mlink] = locdist; + membraneDist[2*mlink+1] = dist; mlink++; } } diff --git a/analysis/Membrane.h b/analysis/Membrane.h index 3d64899a..98c64715 100644 --- a/analysis/Membrane.h +++ b/analysis/Membrane.h @@ -21,6 +21,9 @@ class Membrane { public: + int *membraneLinks; // D3Q19 links that cross membrane + double *membraneDist; // distance to membrane for each linked site + /** * \brief Create a flow adaptor to operate on the LB model * @param ScaLBL - originating data structures @@ -45,9 +48,7 @@ public: * @param membrane - links that form the membrane * @param Np - number of sites in compact layout */ - int Create(DoubleArray &Distance, IntArray &Map, int *neighborList, int *membrane, int Np); - - + int Create(DoubleArray &Distance, IntArray &Map, int *neighborList, int Np); private: }; diff --git a/common/ScaLBL.cpp b/common/ScaLBL.cpp index e3c3baf7..80f39033 100644 --- a/common/ScaLBL.cpp +++ b/common/ScaLBL.cpp @@ -1,5 +1,5 @@ /* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University + Copyright 2013--2022 James E. McClure, Virginia Polytechnic & State University Copyright Equnior ASA This file is part of the Open Porous Media project (OPM). @@ -15,10 +15,8 @@ along with OPM. If not, see . */ #include "common/ScaLBL.h" - #include - ScaLBL_Communicator::ScaLBL_Communicator(std::shared_ptr Dm){ //...................................................................................... Lock=false; // unlock the communicator @@ -324,7 +322,6 @@ ScaLBL_Communicator::ScaLBL_Communicator(std::shared_ptr Dm){ } - ScaLBL_Communicator::~ScaLBL_Communicator() { @@ -988,7 +985,6 @@ int ScaLBL_Communicator::MemoryOptimizedLayoutAA(IntArray &Map, int *neighborLis return(Np); } - void ScaLBL_Communicator::SetupBounceBackList(IntArray &Map, signed char *id, int Np, bool SlippingVelBC) { diff --git a/common/ScaLBL.h b/common/ScaLBL.h index 2be0121e..cb376eb9 100644 --- a/common/ScaLBL.h +++ b/common/ScaLBL.h @@ -699,6 +699,14 @@ public: double GetPerformance(int *NeighborList, double *fq, int Np); int MemoryOptimizedLayoutAA(IntArray &Map, int *neighborList, signed char *id, int Np, int width); + /** + * \brief Create membrane data structure + * - cut lattice links based on distance map + * @param Distance - signed distance to membrane + * @param neighborList - data structure that retains lattice links + * @param Np - number of lattice sites + * @param width - halo width for the model + */ void Barrier(){ ScaLBL_DeviceBarrier(); MPI_COMM_SCALBL.barrier();