electro_kinetic_flows into master

staggered_grid/exec/chargedFluid/test/Makefile.Donev

@@ -0,0 +1,103 @@
+# FBOXLIB_HOME defines the directory in which we will find all the amrex code
+# If you set FBOXLIB_HOME as an environment variable, this line will be ignored
+FBOXLIB_HOME := ../../../../../FBoxLib
+HYDROLIB_HOME := ../../../../HydroGrid
+STAGGERED_HOME := ../../..
+
+# Instructions:
+# NDEBUG=t means optimize, not setting it means debug mode
+# MPI=t means use MPI, otherwise compile for serial execution
+# Important: Choose compilers (gfortran or ifort)
+
+NDEBUG    := t
+MPI       := 
+OMP       :=
+PROF      :=
+COMP      := gfortran
+CCOMP     := gcc
+MKVERBOSE := t
+
+# For Intel compiler (requires ifx+icc as of version 2023)
+#COMP      := Intel
+#CCOMP     := icc
+
+# need this to compile bl_rng.f90 (module to store random number engines and distributions)
+CXX11 = t
+
+include $(FBOXLIB_HOME)/Tools/F_mk/GMakedefs.mak
+
+include $(STAGGERED_HOME)/src_lowMach/GPackage.mak
+VPATH_LOCATIONS += $(STAGGERED_HOME)/src_lowMach
+
+include $(STAGGERED_HOME)/src_charged/GPackage.mak
+VPATH_LOCATIONS += $(STAGGERED_HOME)/src_charged
+
+include $(STAGGERED_HOME)/src_multiSpec/GPackage.mak
+VPATH_LOCATIONS += $(STAGGERED_HOME)/src_multiSpec
+
+include $(STAGGERED_HOME)/src_chemistry/GPackage.mak
+VPATH_LOCATIONS += $(STAGGERED_HOME)/src_chemistry
+
+include $(STAGGERED_HOME)/src_gmres/GPackage.mak
+VPATH_LOCATIONS += $(STAGGERED_HOME)/src_gmres
+
+include $(STAGGERED_HOME)/src_common/GPackage.mak
+VPATH_LOCATIONS += $(STAGGERED_HOME)/src_common
+
+include $(HYDROLIB_HOME)/GPackage.mak
+VPATH_LOCATIONS += $(HYDROLIB_HOME)
+
+include $(FBOXLIB_HOME)/Src/BaseLib/GPackage.mak
+VPATH_LOCATIONS += $(FBOXLIB_HOME)/Src/BaseLib
+
+include $(FBOXLIB_HOME)/Src/MultiGrid/GPackage.mak
+VPATH_LOCATIONS += $(FBOXLIB_HOME)/Src/MultiGrid/
+
+# Because module files are not portable, we would need to recompile
+# LAPACK95 unless we use Intel for everything.
+# Therefore, Donev has disabled LAPACK95
+#LAPACK95_DIR = $(STAGGERED_HOME)/src_multiSpec/LAPACK95
+#F90FLAGS += -I$(LAPACK95_DIR)/lapack95_modules
+
+# Linking lapack and blas, either default ones or Intel ones
+# Intel OneAPI is assumed to be installed, which it is on both CIMS and donevm laptop
+# For Intel, best is to just link this dynamic "use all of MKL"
+# $(MKLROOT)/lib/intel64/libmkl_rt.so
+# If you want static then link in
+# $(MKLROOT)/lib/intel64/libmkl_blas95_lp64.a
+# $(MKLROOT)/lib/intel64/libmkl_lapack95_lp64.a 
+
+ifeq ($(findstring cims.nyu.edu, $(HOSTNAME)), cims.nyu.edu)
+   # Best to use Intel's LAPACK/BLAS on CIMS to avoid old gcc versions in the mix
+   # Make sure to "source /opt/pkg/intel/oneapi/setvars.sh" for ifort (module load intel-oneapi)
+   # This exports variables $(MKLROOT) and $(I_MPI_ROOT)
+   #MKLROOT=/opt/pkg/intel/oneapi/mkl/latest
+else ifeq ($(findstring pop-os, $(HOSTNAME)), pop-os) 
+   # ubuntu system76 pop-os   
+   # Make sure to "source /home/donev/intel/oneapi/setvars.sh" for ifort
+   # MKLROOT=/home/donev/intel/oneapi/mkl/latest
+endif
+
+ifeq ($(findstring gfortran, $(FC)), gfortran)
+   # Since LAPACK/BLAS are compiled with old gcc version
+   # Use instead the Intel MPI library
+   ifdef MKLROOT
+      # Use the single-library dynamic interface of Intel meant for GNU compiler:
+      libraries += $(MKLROOT)/lib/intel64/libmkl_rt.so      
+      MKL_INTERFACE_LAYER=MKL_INTERFACE_LP64+MKL_INTERFACE_GNU
+   else
+      libraries += -llapack -lblas
+   endif   
+else ifeq ($(findstring if, $(FC)), if) # ifx or ifort (Intel)
+   # The Intel compiler mkl flag should do everything on its own, including linking:
+   F90FLAGS += -qmkl
+   MKL_INTERFACE_LAYER=MKL_INTERFACE_LP64
+else
+   libraries += -llapack -lblas     
+endif
+
+main.$(suf).exe: $(objects) 
+	$(LINK.f90) -o main.$(suf).exe $(objects) $(libraries)
+
+include $(FBOXLIB_HOME)/Tools/F_mk/GMakerules.mak
+

staggered_grid/exec/chargedFluid/test/README

@@ -1,3 +1,11 @@
+NEW -- two directories 
+
+electro_osmosis_case1
+electro_osmosis_case2
+
+are examples that run by reading in a pnp profile
+
+
 Supported Problems (./inputs_files/):
 
 ----------------------------------

staggered_grid/exec/chargedFluid/test/eo_post_proc_python_scripts/debye_huckel_check.py

@@ -0,0 +1,107 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+debye_len_case = '1_50'
+Ly = 1.28e-4
+#debye_len_case = '1_10'
+
+if (debye_len_case == '1_50'):
+	E_ext = 3.e9
+	dphi_dn_wall = 794201762623.94604
+elif (debye_len_case == '1_10'): 
+	E_ext = 3.e10
+	dphi_dn_wall = 31767657979.612209
+
+# import deterministic potential and velocity
+
+u_det = np.loadtxt('u_deterministic_' + debye_len_case + '.txt', unpack = True, usecols = [0])
+epot = np.loadtxt('epot_deterministic_' + debye_len_case + '.txt', unpack = True, usecols = [0])
+c1 = np.loadtxt('c1_deterministic_' + debye_len_case + '.txt', unpack = True, usecols = [0])
+c2 = np.loadtxt('c2_deterministic_' + debye_len_case + '.txt', unpack = True, usecols = [0])
+c3 = np.loadtxt('c3_deterministic_' + debye_len_case + '.txt', unpack = True, usecols = [0])
+
+
+#####################
+##  need to calculuate lambda_d based on centerline concentration values: 
+#####################
+mu = 1.05e-2                           # momentum diffusion coefficient
+dielectric_const = 6.91e-19            # dielectric constant
+k_B = 1.3806488e-16
+T = 300.
+
+q1 = 4.2e3
+q2 = -2.72e3 
+
+m1 = 3.82e-23
+m2 = 5.89e-23
+m3 = 3.35e-23
+rho = 1.
+c1_c = c1[len(c1)/2]
+c2_c = c2[len(c2)/2]
+c3_c = c3[len(c3)/2]
+epot_c = epot[len(epot)/2]
+u_c = u_det[len(u_det)/2]
+
+######################
+## check concentration profile against Boltzmann distribution
+######################
+c1_boltz = c1_c*np.exp(-m1*q1/k_B/T*(epot-epot_c))
+c1_approx = c1_c*(1. - m1*q1/k_B/T*(epot-epot_c))
+
+c2_boltz = c2_c*np.exp(-m2*q2/k_B/T*(epot-epot_c))
+c2_approx = c2_c*(1. - m2*q2/k_B/T*(epot-epot_c))
+
+
+#####################
+##  check relationship: 
+##
+##  u = \epsilon E_x/ \mu (\phi - \phi_{wall})
+## 
+#####################
+shift = epot_c - u_c*mu/dielectric_const/E_ext
+epot_wall = epot[0]
+u_theory = dielectric_const*E_ext/mu * (epot - shift)
+
+######################
+## compute debye length based on concentrations at channel center
+######################
+lambda_d = np.sqrt(dielectric_const*k_B*T/(rho*(c1_c*m1*q1**2 + c2_c*m2*q2**2)))
+print 'debye len is: ', lambda_d
+print 'ratio Ly/lambda_d: ', Ly/lambda_d
+
+
+
+
+######################
+## compare centerline velocity with D-H theory
+######################
+print 'Debye Huckel slip velocity: ', lambda_d*dielectric_const*dphi_dn_wall*E_ext/mu
+print 'Deterministic calculation vel: ', u_c
+
+
+
+plt.figure(1)
+plt.plot(u_det, 'b-o', label = 'simulation $u$')
+plt.plot(u_theory, 'r-o', label = '$\epsilon E_x (\phi-\phi_w)/\mu$')
+plt.title('Velocity')
+plt.legend(loc = 'best')
+
+plt.figure(5)
+plt.plot(u_det-u_theory, 'g-o')
+plt.title('Difference between velocities')
+
+plt.figure(2)
+plt.plot(c1, 'g-o', label = 'simulation')
+plt.plot(c1_boltz, 'r-o', label = 'theory') 
+plt.plot(c1_approx, 'b-o', label = 'approximation') 
+plt.title('Concentration 1')
+plt.legend(loc='best')
+
+plt.figure(3)
+plt.plot(c2, 'g-o', label = 'simulation')
+plt.plot(c2_boltz, 'r-o', label = 'theory') 
+plt.plot(c2_approx, 'b-o', label = 'approximation') 
+plt.title('Concentration 2')
+plt.legend(loc='best')
+
+plt.show()

staggered_grid/exec/chargedFluid/test/eo_post_proc_python_scripts/plot_1d_profs.py

@@ -0,0 +1,76 @@
+import os
+import urllib2
+import numpy as np
+import matplotlib.font_manager as font_manager
+import matplotlib.pyplot as plt
+import pylab
+
+
+##==========================================
+## Get time averaged profiles
+##
+
+#debye_len_case = '1_50'
+Ly = 1.28e-4
+debye_len_case = '1_10'
+
+if (debye_len_case == '1_50'): 
+	y_vals, u, v, w, c1, c2, c3, epot = np.loadtxt('hstat00064000_debye_len_1_50', unpack = True, usecols = [0,1,2,3,5,6,7,9])
+elif (debye_len_case == '1_10'): 
+	y_vals, u, v, w, c1, c2, c3, epot = np.loadtxt('hstat00056500_debye_len_1_10', unpack = True, usecols = [0,1,2,3,5,6,7,9])
+
+u_det = np.loadtxt('u_deterministic_' + debye_len_case + '.txt', unpack = True, usecols = [0])
+c1_det = np.loadtxt('c1_deterministic_' + debye_len_case + '.txt', unpack = True, usecols = [0])
+c2_det = np.loadtxt('c2_deterministic_' + debye_len_case + '.txt', unpack = True, usecols = [0])
+c3_det = np.loadtxt('c3_deterministic_' + debye_len_case + '.txt', unpack = True, usecols = [0])
+epot_det = np.loadtxt('epot_deterministic_' + debye_len_case + '.txt', unpack = True, usecols = [0])
+
+
+##==========================================
+## Plot results 
+##
+
+plt.figure(1) 
+plt.plot(y_vals, u, 'b-o', label = '3d averaged')
+plt.plot(y_vals, u_det, 'r-o', label = '2d')
+plt.title('$u$ v. wall-normal $y$')
+plt.legend(loc='best')
+
+plt.figure(7) 
+plt.plot(y_vals, v, 'b-o', label = '3d averaged')
+plt.title('$v$ v. wall-normal $y$')
+plt.legend(loc='best')
+
+plt.figure(9) 
+plt.plot(y_vals, w, 'b-o', label = '3d averaged')
+plt.title('$w$ v. wall-normal $y$')
+plt.legend(loc='best')
+
+plt.figure(2) 
+plt.plot(y_vals, c1, 'b-o', label = '3d averaged')
+plt.plot(y_vals, c1_det, 'r-o', label = '2d')
+plt.title('$c_1$ v. wall-normal $y$')
+plt.legend(loc='best')
+
+plt.figure(3) 
+plt.plot(y_vals, c2, 'b-o', label = '3d averaged')
+plt.plot(y_vals, c2_det, 'r-o', label = '2d')
+plt.title('$c_2$ v. wall-normal $y$')
+plt.legend(loc='best')
+
+plt.figure(4) 
+plt.plot(y_vals, c3, 'b-o', label = '3d averaged')
+plt.plot(y_vals, c3_det, 'r-o', label = '2d')
+plt.title('$c_3$ v. wall-normal $y$')
+plt.legend(loc='best')
+
+plt.figure(5) 
+plt.plot(y_vals, epot, 'b-o', label = '3d averaged')
+plt.plot(y_vals, epot_det, 'r-o', label = '2d')
+plt.title('$\phi$ v. wall-normal $y$')
+plt.legend(loc='best')
+
+### show ###
+plt.show()
+
+