add monitoring debug codes

src_lbm/shan-chen/GNUmakefile

@@ -1,5 +1,5 @@
-AMREX_HOME ?= ../../../amrex
-FHDEX_HOME ?= ../../../FHDeX
+AMREX_HOME ?= ../../../Code/amrex
+FHDEX_HOME ?= ../../../Code/FHDeX
 
 DEFINES += -DMAX_SPECIES=2
 
@@ -12,7 +12,7 @@ COMP    = gnu
 USE_FFT   = TRUE
 USE_MPI   = FALSE
 USE_OMP   = FALSE
-USE_CUDA  = TRUE
+USE_CUDA  = FALSE
 USE_HIP   = FALSE
 
 include $(AMREX_HOME)/Tools/GNUMake/Make.defs

src_lbm/shan-chen/LBM_Debug.H

@@ -0,0 +1,41 @@
+#ifndef LBM_DEBUG_H_
+#define LBM_DEBUG_H_
+
+#include "LBM_shan-chen.H"
+
+void PrintMFComponent(MultiFab& mfab, int compIdx=0, int ngrow=-1){
+    int ncomp = mfab.nComp();
+    auto const & mfab_ptr = mfab.arrays();
+    if(compIdx >= ncomp) {
+        Print() << "Error: compIdx >= ncomp (" << compIdx << " >= " << ncomp << ") in function [PrintMFComponent] \n";
+        return;
+    }
+    IntVect halo;
+    if(ngrow<0){
+        halo = IntVect(mfab.nGrow());
+    }else{
+        halo = IntVect(ngrow);
+    }
+    // loop over all boxes in the MultiFab
+    for (MFIter mfi(mfab); mfi.isValid(); ++mfi){
+        const Box& bx = grow(mfi.validbox(), halo);
+        const Array4<Real>& fab_arr = mfab.array(mfi);
+        const auto lo = lbound(bx);
+        const auto hi = ubound(bx);
+
+        Print() << "====== Box " << mfi.index() << " component " << compIdx << " ======\n";
+        for (int k = lo.z; k <= hi.z; ++k) {
+            for (int j = lo.y; j <= hi.y; ++j) {
+                for (int i = lo.x; i <= hi.x; ++i) {
+                    Print() << "(" << i << "," << j << "," << k << ") = " << fab_arr(i,j,k,compIdx) << "  ";
+                }
+                Print() << "\n";
+            }
+            Print() << "\n";
+        }
+        Print() << "===========================\n";
+    }
+}
+
+
+#endif

src_lbm/shan-chen/LBM_shan-chen.H

@@ -2,19 +2,23 @@
 #define LBM_H_
 
 #include "LBM_d3q19.H"
+#include "LBM_Debug.H"
+
+const IntVect print_grid(0, 0, 0); // grid to print out
 
 // default initial condition
 std::string init_cond = "mixture";
+const bool use_SC_pseudo = true; // false: use density field directly; true: use SC pseudopotentials;
 
-AMREX_GPU_MANAGED Real G = 1.5;
-AMREX_GPU_MANAGED Real rho_hi = 2.3; // values reported by Belardinelli et al. [Phys. Rev. E 91, 023313 (2015)]
-AMREX_GPU_MANAGED Real rho_lo = 0.06;
+AMREX_GPU_MANAGED Real G = 0.24; //1.5;
+AMREX_GPU_MANAGED Real rho_hi = 0.7; //3.0; // values reported by Belardinelli et al. [Phys. Rev. E 91, 023313 (2015)]
+AMREX_GPU_MANAGED Real rho_lo = 0.0;
 AMREX_GPU_MANAGED Real sigma = 0.174; // interfacial tension
-AMREX_GPU_MANAGED Real alpha = 6; // interface width
+AMREX_GPU_MANAGED Real alpha = 1.; // interface width
 
-AMREX_GPU_MANAGED Real temperature = 1e-5;
+AMREX_GPU_MANAGED Real temperature = 0.; //1e-5;
 
-AMREX_GPU_MANAGED Real init_frac = 0.5; // fractional size of droplet or stripe
+AMREX_GPU_MANAGED Real init_frac = 0.4; // fractional size of droplet or stripe
 
 AMREX_GPU_CONSTANT Real const lambda_d = 1.0;
 AMREX_GPU_CONSTANT Real const lambda_b = 1.0;
@@ -45,8 +49,9 @@ RealVect gradient(int x, int y, int z, const Array4<Real>& field, int icomp) {
     int xp = x + c[i][0];
     int yp = y + c[i][1];
     int zp = z + c[i][2];
+    Real SC_pseudo_neighbor = use_SC_pseudo ? ::rho_hi*(1.-exp(-field(xp,yp,zp,icomp)/::rho_hi)) : field(xp,yp,zp,icomp);
     for (int dir=0; dir<3; dir++) {
-      gradient[dir] += w[i]/cs2*field(xp,yp,zp,icomp)*c[i][dir];
+      gradient[dir] += w[i]/cs2*SC_pseudo_neighbor*c[i][dir];
     }
   }
   return gradient;
@@ -55,17 +60,19 @@ RealVect gradient(int x, int y, int z, const Array4<Real>& field, int icomp) {
 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
 Real laplacian(int x, int y, int z, const Array4<Real>& field, int icomp) {
   Real laplacian = 0.0;
+  Real SC_pseudo_center = use_SC_pseudo ? ::rho_hi*(1.-exp(-field(x,y,z,icomp)/::rho_hi)) : field(x,y,z,icomp);
   for (int i=0; i<nvel; i++) {
     int xp = x + c[i][0];
     int yp = y + c[i][1];
     int zp = z + c[i][2];
-    laplacian += 2.*w[i]/cs2*(field(xp,yp,zp,icomp)-field(x,y,z,icomp));
+    Real SC_pseudo_neighbor = use_SC_pseudo ? ::rho_hi*(1.-exp(-field(xp,yp,zp,icomp)/::rho_hi)) : field(xp,yp,zp,icomp);
+    laplacian += 2.*w[i]/cs2*(SC_pseudo_neighbor-SC_pseudo_center);
   }
   return laplacian;
 }
 
 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
-Array<Array1D<Real,0,nvel>,2> thermal_noise(const Array<Real,2>& rho, const RandomEngine& engine) {
+Array<Array1D<Real,0,nvel>,2> thermal_noise(const Array<Real,2>& rho, const RandomEngine& engine, bool if_print_detail=false) {
   Array<Array1D<Real,0,nvel>,2> xi;
   Array1D<Real,0,nvel> & xif = xi[0];
   Array1D<Real,0,nvel> & xig = xi[1];
@@ -76,23 +83,40 @@ Array<Array1D<Real,0,nvel>,2> thermal_noise(const Array<Real,2>& rho, const Rand
   // fluxes are anticorrelated (total momentum is conserved)
   for (int i=1; i<=AMREX_SPACEDIM; ++i ) {
     const Real pref = (2.-lambda_d)*lambda_d*temperature*rho[0]*rho[1]/(rho[0]+rho[1]);
-    xif(i) = RandomNormal(0., std::sqrt(pref), engine);
+    /*if(if_print_detail) {
+      printf("pref for flux mode %d = %f\n and sqrt(pref) = %f\n", i, pref, std::sqrt(pref));
+    }*/
+    xif(i) = RandomNormal(0., std::sqrt(std::abs(pref)), engine);
     xig(i) = - xif(i);
   }
 
   // remaining modes
   for (int i=AMREX_SPACEDIM+1; i<nvel; ++i) {
     const Real pref = (2.-lambda_f[i])*lambda_f[i]*b[i]*temperature/cs2*rho[0];
     const Real preg = (2.-lambda_g[i])*lambda_g[i]*b[i]*temperature/cs2*rho[1];
-    xif(i) = RandomNormal(0., std::sqrt(pref), engine);
-    xig(i) = RandomNormal(0., std::sqrt(preg), engine);
+    xif(i) = RandomNormal(0., std::sqrt(std::abs(pref)), engine);
+    xig(i) = RandomNormal(0., std::sqrt(std::abs(preg)), engine);
+  }
+
+  if(if_print_detail) {
+    printf("xi_f in thermal_noise() = ");
+    for (int i=0; i<nvel; ++i) {
+      printf("%f\t", xif(i));
+    }
+    printf("\n");
+    printf("xi_g in thermal_noise() = ");
+    for (int i=0; i<nvel; ++i) {
+      printf("%f\t", xig(i));
+    }
+    printf("\n");
   }
 
   return xi;
 }
 
 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
-Array<RealVect,2> transform_velocity(Array<Real,2> const& rho, Array<RealVect,2> const& v, Array<RealVect,2> const& a, Array<Array1D<Real,0,nvel>,2> const& xia) {
+Array<RealVect,2> transform_velocity(Array<Real,2> const& rho, Array<RealVect,2> const& v, Array<RealVect,2> const& a, Array<Array1D<Real,0,nvel>,2> const& xia,
+  bool if_print_detail=false) {
   Array<RealVect, 2> vr;
   RealVect xi;
 
@@ -101,7 +125,13 @@ Array<RealVect,2> transform_velocity(Array<Real,2> const& rho, Array<RealVect,2>
   for (int i=0; i<2; ++i) {
     int j = (i+1)%2;
     for (int k=0; k<AMREX_SPACEDIM; ++k) xi[k] = xia[i](1+k);
-    vr[i] = v[i] + 0.5*(a[i] - lambda_d*rho[j]/(rho[i]+rho[j])*(v[i]-v[j]+0.5*(a[i]-a[j])) + xi/rho[i]);
+    vr[i] = v[i] + 0.5*(a[i] - lambda_d*rho[j]/(rho[i]+rho[j])*(v[i]-v[j]+0.5*(a[i]-a[j])) + ((std::abs(rho[i]) > FLT_EPSILON) ? xi/rho[i] : RealVect(0.)));
+  }
+  if(if_print_detail) {
+    printf("xi in transform_velocity() = %f\t%f\n", xi[0], xi[1]);
+    printf("density in transform_velocity() = %f\t%f\n", rho[0], rho[1]);
+    printf("velocity in transform_velocity() = vf=(%f,%f,%f),\tvg=(%f,%f,%f)\n", vr[0][0], vr[0][1], vr[0][2], vr[1][0], vr[1][1], vr[1][2]);
+    printf("acceleration in transform_velocity() = af=(%f,%f,%f),\tag=(%f,%f,%f)\n", a[0][0], a[0][1], a[0][2], a[1][0], a[1][1], a[1][2]);
   }
 
   return vr;
@@ -160,7 +190,7 @@ void collide(int x, int y, int z,
   rho_ref[0] = ref(x,y,z,0);
   rho_ref[1] = ref(x,y,z,1);
 
-  xi = thermal_noise(rho_ref, engine);
+  // xi = thermal_noise(rho_ref, engine);
   Array1D<Real,0,nvel> const& xif = xi[0];
   Array1D<Real,0,nvel> const& xig = xi[1];
 
@@ -170,21 +200,41 @@ void collide(int x, int y, int z,
   rho[0] = mf(0);
   rho[1] = mg(0);
 
+  if(x==print_grid[0] && y==print_grid[1] && z==print_grid[2]){
+    xi = thermal_noise(rho, engine /*, true*/);  // use current densities for noise calculation
+  }else{
+    xi = thermal_noise(rho, engine);
+  }
+
   v[0] = (mf(0) > FLT_EPSILON) ? RealVect(mf(1), mf(2), mf(3))/mf(0) : RealVect(0.);
   v[1] = (mg(0) > FLT_EPSILON) ? RealVect(mg(1), mg(2), mg(3))/mg(0) : RealVect(0.);
 
-  a[0] = -G*cs2*gradient(x,y,z,h,1);
-  a[1] = -G*cs2*gradient(x,y,z,h,0);
+  Real SC_pseudo_rho = use_SC_pseudo ? ::rho_hi*(1.-exp(-rho[0]/::rho_hi)) : rho[0];
+  Real SC_pseudo_phi = use_SC_pseudo ? ::rho_hi*(1.-exp(-rho[1]/::rho_hi)) : rho[1];
+  a[0] = (std::abs(rho[0]) > FLT_EPSILON)? -G*cs2*SC_pseudo_rho*gradient(x,y,z,h,1)/rho[0] : RealVect(0.);
+  a[1] = (std::abs(rho[1]) > FLT_EPSILON)? -G*cs2*SC_pseudo_phi*gradient(x,y,z,h,0)/rho[1] : RealVect(0.);
 
-  vr = transform_velocity(rho, v, a, xi);
+  if(x==print_grid[0] && y==print_grid[1] && z==print_grid[2]){
+    vr = transform_velocity(rho, v, a, xi /*, true*/);
+  }else{
+    vr = transform_velocity(rho, v, a, xi);
+  }
 
   v_b = (rho[0]*v[0]+rho[1]*v[1])/(rho[0]+rho[1]);
   a_b = (rho[0]*a[0]+rho[1]*a[1])/(rho[0]+rho[1]);
 
   mfEq = equilibrium_moments(rho[0], v_b+0.5*a_b);
   mgEq = equilibrium_moments(rho[1], v_b+0.5*a_b);
 
-  phi = force_moments(rho, vr, a);
+  /*if(x==print_grid[0] && y==print_grid[1] && z==print_grid[2]){
+    //printf("fluid f,g velocity in collide() = \n");
+    printf("xi before force_moments = %f\t%f\t%f\n", xi[0](0), xi[0](1), xi[0](2));
+    phi = force_moments(rho, vr, a);
+    printf("\n");
+  }else{*/
+    phi = force_moments(rho, vr, a);
+  //}
+
   Array1D<Real,0,nvel> const& phif = phi[0];
   Array1D<Real,0,nvel> const& phig = phi[1];
 
@@ -310,7 +360,7 @@ void hydrovars(int x, int y, int z,
   Array<Real,2> rho;
   Array<RealVect,2> v, a, vr;
 
-  rho[0] = mf(0);
+  rho[0] = mf(0); 
   rho[1] = mg(0);
 
   v[0] = (mf(0) > FLT_EPSILON) ? RealVect(mf(1), mf(2), mf(3))/mf(0) : RealVect(0.);
@@ -450,9 +500,11 @@ void LBM_init_droplet(const Real dfrac,
 void LBM_init(const Geometry& geom, MultiFab& mf, MultiFab& mg, MultiFab& hydrovs, MultiFab& refstate) {
   if (init_cond == "mixture") LBM_init_mixture(geom, mf, mg, hydrovs, refstate);
   else if (init_cond == "stripe") LBM_init_stripe(init_frac, geom, mf, mg, hydrovs, refstate);
-  else if (init_cond == "droplet") LBM_init_droplet(init_frac, geom, mf, mg, hydrovs, refstate);
-  else {
-    Print() << "ERROR: Unkown init_cond";
+  else if (init_cond == "droplet") {
+    LBM_init_droplet(init_frac, geom, mf, mg, hydrovs, refstate);
+    Print() << "Initialized droplet with frac radius = " << init_frac/2.0 << "\n";
+  } else {
+    Print() << "ERROR: Unknown init_cond";
     ParallelDescriptor::Abort(-1, false);
   }
 }

src_lbm/shan-chen/Make.package

@@ -1,3 +1,7 @@
 CEXE_headers += LBM_d3q19.H
 CEXE_headers += LBM_binary.H
+CEXE_headers += LBM_shan-chen.H
 CEXE_headers += LBM_fluctuations.H
+CEXE_headers += LBM_Debug.H
+
+CEXE_sources += main_driver.cpp

src_lbm/shan-chen/Viewer.ipynb

@@ -0,0 +1,105 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "30285ee4",
+   "metadata": {},
+   "source": [
+    "### Generate droplet trajectories\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b2836b93",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import yt\n",
+    "import os\n",
+    "import numpy as np\n",
+    "import glob\n",
+    "import re\n",
+    "\n",
+    "# Parameters\n",
+    "alpha0 = 1.5\n",
+    "temperature = 0.0\n",
+    "nx, ny, nz = 32, 32, 32\n",
+    "radius = 0.2\n",
+    "rho_lo = 0.0\n",
+    "rho_hi = 3.0\n",
+    "step1 = 360000\n",
+    "step2 = 600000\n",
+    "plot_int = 1000\n",
+    "\n",
+    "# Construct file paths\n",
+    "data_dir = \"./data_droplet_alpha0_{:.2f}_r{:.2f}_size{:d}-{:d}-{:d}/\".format(alpha0, radius, nx, ny, nz)\n",
+    "\n",
+    "IMG_folder = data_dir + \"/IMG_data_shshan_alpha0_{:.2f}_xi_{:.1e}_size{:d}-{:d}-{:d}/\".format(alpha0, temperature, nx, ny, nz)\n",
+    "\n",
+    "# Create image directory if it doesn't exist\n",
+    "os.makedirs(IMG_folder, exist_ok=True)\n",
+    "all_plt_files = sorted(glob.glob(os.path.join(data_dir, \"plt*\")))\n",
+    "\n",
+    "# Extract step numbers from filenames (assuming format like plt00010, plt12345, etc.)\n",
+    "def extract_step(filepath):\n",
+    "    basename = os.path.basename(filepath) # Get the file/folder name from the path\n",
+    "    match = re.search(r'plt(\\d+)', basename)  # It finds single digits with \\d and sequences of digits with \\d+.\n",
+    "    if match:\n",
+    "        return int(match.group(1)) # Extract the digits in filename and convert to integer\n",
+    "    else:\n",
+    "        return -1  # invalid\n",
+    "# Load all plot files (e.g., plt0000, plt0001, ...)\n",
+    "ts = yt.load(data_dir + \"plt*\")\n",
+    "\n",
+    "# Filter files within [step1, step2] and divisible by plot_int (or aligned to step1 + n*plot_int)\n",
+    "selected_files = []\n",
+    "for f in all_plt_files:\n",
+    "    step = extract_step(f)\n",
+    "    if step == -1:\n",
+    "        continue\n",
+    "    if step1 <= step <= step2 and (step - step1) % plot_int == 0:\n",
+    "        selected_files.append((step, f))\n",
+    "\n",
+    "# Sort by step number\n",
+    "selected_files.sort(key=lambda x: x[0])\n",
+    "\n",
+    "# Process each selected file\n",
+    "ac = 0\n",
+    "for step, filepath in selected_files:\n",
+    "    fr = yt.load(filepath)\n",
+    "    slc = yt.SlicePlot(fr, \"z\", (\"boxlib\", \"rhoA\"))\n",
+    "    slc.set_log((\"boxlib\", \"rhoA\"), False)\n",
+    "    # Optional: set color limits\n",
+    "    # slc.set_zlim((\"boxlib\", \"rho\"), rho_lo - 0.3, rho_hi + 0.5)\n",
+    "    \n",
+    "    filename = os.path.join(IMG_folder, f\"fr{ac:07d}.png\")\n",
+    "    slc.save(filename)\n",
+    "    ac += 1\n",
+    "\n",
+    "print(f\"Saved {ac} frames from step {step1} to {step2} (interval={plot_int})\")"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "amrex_env",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.13.5"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}