BinWang0213/mpi_navierstokes.py

## readme.md

      
    Raw
  

              readme.md
            
          
    ngsolve 2021.5.5 install on LSU QueenBee-3 Cluster with MKL


Follow the guide @ https://github.com/mcovalt/FEniCS-MKL/blob/master/Dockerfile
HPC command https://info.hpc.sussex.ac.uk/hpc-guide/resource.html
guide2 https://wiki.math.uwaterloo.ca/fluidswiki/index.php?title=FEniCS

Enter a node ssh mike012
Login:
ssh -X bwang31@qbc.loni.org
srun -t 5:00:00 -n48 -N1 -A loni_mp2020 -p bigmem --pty /bin/bash
srun -t 5:00:00 -n48 -N1 -A loni_mp2020 -p checkpt --pty /bin/bash


MPI and compiling env setup

#Using node for faster compile
srun -t 4:00:00 -n24 -N1 -A loni_mp2020 -p workq --pty /bin/bash

export NGS_ROOT=/project/karsten/bwang31/ngsolve
mkdir -p $NGS_ROOT && cd $NGS_ROOT

module purge
module load python/3.7.6
module load cmake/3.16.2/intel-19.0.5
module load gcc/9.3.0

export CFLAGS='-O3'
export CXXFLAGS='-O3'
export FCFLAGS='-O3'
export F77FLAGS='-O3'
export F90FLAGS='-O3'
export MPICCFLAGS='-O3'
export MPICXXFLAGS='-O3'
export MPIF77FLAGS='-O3'
export MPIF90FLAGS='-O3'
export MKL_NUM_THREADS=1
export OMP_NUM_THREADS=1
Build openMPI

cd $NGS_ROOT
wget -P download "https://download.open-mpi.org/release/open-mpi/v3.1/openmpi-3.1.6.tar.gz"

mkdir -p deps && cd deps
tar -zxf ../download/openmpi-3.1.6.tar.gz

cd $NGS_ROOT/deps/openmpi-3.1.6
#without-verbs is used to disable " OpenFabrics device error "
#ref: https://www.open-mpi.org/faq/?category=openfabrics
./configure --prefix=$NGS_ROOT/deps/openmpi-3.1.6/openmpi-lib
make all install

export OPENMPI_PATH=$NGS_ROOT/deps/openmpi-3.1.6/openmpi-lib
export PATH=$PATH:$OPENMPI_PATH/bin
export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$OPENMPI_PATH/lib

Python environment setup

export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/home/packages/python/3.7.6/lib/

#small dependency for ngsolve
pip install --user mpi4py
pip install petsc petsc4py #this is only required for ngsolve-petsc module

Build ngsolve

export BASEDIR=$NGS_ROOT
cd $BASEDIR && git clone https://github.com/NGSolve/ngsolve.git ngsolve-src
cd $BASEDIR/ngsolve-src && git submodule update --init --recursive
mkdir $BASEDIR/ngsolve-build
mkdir $BASEDIR/ngsolve-install

#We need to define this first to avoid error while building
export NETGENDIR="${BASEDIR}/ngsolve-install/bin"
export PATH=$NETGENDIR:$PATH
export PYTHONPATH=$NETGENDIR/../`python3 -c "from distutils.sysconfig import get_python_lib; print(get_python_lib(1,0,''))"`:$PYTHONPATH

cd $BASEDIR/ngsolve-build
cmake -DCMAKE_INSTALL_PREFIX=${BASEDIR}/ngsolve-install ${BASEDIR}/ngsolve-src \
  -DCMAKE_BUILD_TYPE=RELEASE \
  -DBUILD_STUB_FILES=OFF \
  -DPYTHON_EXECUTABLE:FILEPATH=$(which python) \
  -DPYTHON_INCLUDE_DIR=$(python -c "from distutils.sysconfig import get_python_inc; print(get_python_inc())")  \
  -DPYTHON_LIBRARY=$(python -c "import distutils.sysconfig as sysconfig; print(sysconfig.get_config_var('LIBDIR'))") \
  -DUSE_GUI=OFF \
  -DUSE_UMFPACK=ON \
  -DUSE_MUMPS=ON \
  -DUSE_HYPRE=OFF \
  -DUSE_MPI=ON \
  -DMKL_SDL=OFF \
  -DUSE_MKL=ON \
  -DLAPACK_LIBRARIES=${MKLROOT}/lib/intel64 \
  -DCMAKE_CXX_FLAGS="-O3 -std=c++17" \
  -DCMAKE_C_FLAGS="-O3"

make # parallel build may fail
make install
ngsolve environment

/project/karsten/bwang31/ngsolve/bashrc.sh
export NGS_ROOT=/project/karsten/bwang31/ngsolve

module purge
module load python/3.7.6
module load cmake/3.16.2/intel-19.0.5
module load gcc/9.3.0

export CFLAGS='-O3'
export CXXFLAGS='-O3'
export FCFLAGS='-O3'
export F77FLAGS='-O3'
export F90FLAGS='-O3'
export MPICCFLAGS='-O3'
export MPICXXFLAGS='-O3'
export MPIF77FLAGS='-O3'
export MPIF90FLAGS='-O3'
export MKL_NUM_THREADS=1
export OMP_NUM_THREADS=1

export OPENMPI_PATH=$NGS_ROOT/deps/openmpi-3.1.6/openmpi-lib
export PATH=$PATH:$OPENMPI_PATH/bin
export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$OPENMPI_PATH/lib

export BASEDIR=$NGS_ROOT
export NETGENDIR="${BASEDIR}/ngsolve-install/bin"
export PATH=$NETGENDIR:$PATH
export PYTHONPATH=$NETGENDIR/../`python3 -c "from distutils.sysconfig import get_python_lib; print(get_python_lib(1,0,''))"`:$PYTHONPATH

cd /work/bwang31/ngsolve

Testing ngsolve (we need to run ngsolve with mpirun even in a serial mode)

source /project/karsten/bwang31/ngsolve/bashrc.sh

cd /work/bwang31/ngsolve
cp ${BASEDIR}/ngsolve-install/share/ngsolve/py_tutorials/mpi/mpi_navierstokes.py .
# Copy the following modified mpi_naiversoktes

# Run simulation
#--mca btl '^openib' is used to disable openfabrics error
mpirun -np 8 ngspy mpi_navierstokes.py --mca btl '^openib'

#Union VTK into a single time-series XDMF

#Must using the following fname format
#file_name = f"{output_path}/vtkout_time{t:.3f}_proc{rank} "
python unionSolution.py --dir navierstokes_output

#View xdmf and h5 using paraview


## mpi_navierstokes.py
import netgen.meshing
from netgen.geom2d import SplineGeometry


from ngsolve import *

comm = mpi_world
rank = comm.rank
np = comm.size

do_vtk = True

# viscosity
nu = 0.001

# timestepping parameters
tau = 0.001
tend = 3.0

# mesh = Mesh("cylinder.vol")
geo = SplineGeometry()
geo.AddRectangle( (0, 0), (2, 0.41), bcs = ("wall", "outlet", "wall", "inlet"))
geo.AddCircle ( (0.2, 0.2), r=0.05, leftdomain=0, rightdomain=1, bc="cyl")
if rank==0:
    ngmesh = geo.GenerateMesh(maxh=0.08)
    ngmesh.Distribute(comm)
else:
    ngmesh = netgen.meshing.Mesh.Receive(comm)
    ngmesh.SetGeometry(geo)

mesh = Mesh(ngmesh)

print(f"Rank{rank} NumEles={mesh.ne}")

#does not work with mpi yet...
#mesh.Curve(3)

V = H1(mesh,order=3, dirichlet="wall|cyl|inlet")
Q = H1(mesh,order=2)

X = FESpace([V,V,Q])

ux,uy,p = X.TrialFunction()
vx,vy,q = X.TestFunction()

div_u = grad(ux)[0]+grad(uy)[1]
div_v = grad(vx)[0]+grad(vy)[1]

stokes = nu*grad(ux)*grad(vx)+nu*grad(uy)*grad(vy)+div_u*q+div_v*p - 1e-10*p*q
a = BilinearForm(X)
a += SymbolicBFI(stokes)
a.Assemble()

# nothing here ...
f = LinearForm(X)
f.Assemble()

# gridfunction for the solution
gfu = GridFunction(X)

# parabolic inflow at bc=1:
uin = 1.5*4*y*(0.41-y)/(0.41*0.41)
gfu.components[0].Set(uin, definedon=mesh.Boundaries("inlet"))

velocity = CoefficientFunction(gfu.components[0:2])


# solve Stokes problem for initial conditions:
#inv_stokes = a.mat.Inverse(X.FreeDofs(), inverse="mumps")
inv_stokes = a.mat.Inverse(X.FreeDofs(), inverse="masterinverse")
res = f.vec.CreateVector()
res.data = f.vec - a.mat*gfu.vec
gfu.vec.data += inv_stokes * res

# matrix for implicit Euler
mstar = BilinearForm(X)
mstar += SymbolicBFI(ux*vx+uy*vy + tau*stokes)
mstar.Assemble()

# inv = mstar.mat.Inverse(X.FreeDofs(), inverse="masterinverse")
inv = mstar.mat.Inverse(X.FreeDofs(), inverse="masterinverse")

# the non-linear term
conv = BilinearForm(X, nonassemble = True)
conv += SymbolicBFI( CoefficientFunction( (ux,uy) ) * (grad(ux)*vx+grad(uy)*vy) )

t = 0
vtk_interval = int(0.05/tau);

import os
#output_path = os.path.dirname(os.path.realpath(__file__)) + "/navierstokes_output"
output_path = os.path.dirname(os.path.realpath(__file__)) + "/navierstokes_output"
if rank==0 and not os.path.exists(output_path):
    os.mkdir(output_path)
comm.Barrier() #wait until master has created the directory!!

t=0.0

if rank>0 and do_vtk:
    file_name = f"{output_path}/vtkout_time{t:.3f}_proc{rank} "
    print("rank "+str(rank)+", output to "+file_name)
    vtk = VTKOutput(ma=mesh,coefs=[velocity],names=["u"],filename=file_name,subdivision=1)
    vtk.Do()

count = 1
# implicit Euler/explicit Euler splitting method:
with TaskManager():
    while t < tend:
        if rank==0:
            print ("t=", t)

        conv.Apply (gfu.vec, res)
        res.data += a.mat*gfu.vec
        gfu.vec.data -= tau * inv * res

        if rank>0 and count%vtk_interval==0 and do_vtk:
            file_name = f"{output_path}/vtkout_time{t:.3f}_proc{rank} "
            print("rank "+str(rank)+", output to "+file_name)
            vtk = VTKOutput(ma=mesh,coefs=[velocity],names=["u"],filename=file_name,subdivision=1)
            vtk.Do()
        count = count+1

        t = t + tau
        comm.Barrier()

## unionSolution.py
import numpy as np
from numpy.core.shape_base import block
import glob, os
import re
from collections import defaultdict

import meshio
import pyvista
import pyvista as pv

import argparse

parser = argparse.ArgumentParser(description='This is a module to union parallel vtks for ngsolve.')
parser.add_argument("--dir", default=".", help="Number of processor")
args = parser.parse_args()
print("Input Parameters=",args)


try:
    from vtkmodules.vtkCommonCore import vtkVersion
    VTK9 = vtkVersion().GetVTKMajorVersion() >= 9
except ImportError:  # pragma: no cover
    VTK9 = False

def pv2meshio(mesh):
    #save_meshio() in pyvista/utilities/fileio.py
    from meshio.vtk._vtk import vtk_to_meshio_type

    # Cast to pyvista.UnstructuredGrid
    if not isinstance(mesh, pv.UnstructuredGrid):
        mesh = mesh.cast_to_unstructured_grid()

    # Copy useful arrays to avoid repeated calls to properties
    vtk_offset = mesh.offset
    vtk_cells = mesh.cells
    vtk_cell_type = mesh.celltypes

    # Check that meshio supports all cell types in input mesh
    pixel_voxel = {8, 11}       # Handle pixels and voxels
    for cell_type in np.unique(vtk_cell_type):
        if cell_type not in vtk_to_meshio_type.keys() and cell_type not in pixel_voxel:
            raise TypeError(f"meshio does not support VTK type {cell_type}.")

    # Get cells
    cells = []
    c = 0
    for offset, cell_type in zip(vtk_offset, vtk_cell_type):
        numnodes = vtk_cells[offset+c]
        if VTK9:  # must offset by cell count
            cell = vtk_cells[offset+1+c:offset+1+c+numnodes]
            c += 1
        else:
            cell = vtk_cells[offset+1:offset+1+numnodes]
        cell = (
            cell if cell_type not in pixel_voxel
            else cell[[0, 1, 3, 2]] if cell_type == 8
            else cell[[0, 1, 3, 2, 4, 5, 7, 6]]
        )
        cell_type = cell_type if cell_type not in pixel_voxel else cell_type+1
        cell_type = (
            vtk_to_meshio_type[cell_type] if cell_type != 7
            else f"polygon{numnodes}"
        )

        if len(cells) > 0 and cells[-1][0] == cell_type:
            cells[-1][1].append(cell)
        else:
            cells.append((cell_type, [cell]))

    for k, c in enumerate(cells):
        cells[k] = (c[0], np.array(c[1]))


    return  meshio.Mesh(np.array(mesh.points),cells)

#Search vtk files in a given path
os.chdir(args.dir)
vtks=[f for f in glob.glob("*.vtk")]
vtks.sort(key=str.lower)
#print('VtkFiles=',vtks)

#Parse filename and group vtk files based on time
#ngsolve must using the following format
#file_name = f"{output_path}/vtkout_time{t:.3f}_proc{rank} "

vtks_group=defaultdict(list)
for fname in vtks:
    #https://stackoverflow.com/questions/4289331/how-to-extract-numbers-from-a-string-in-python
    rr = re.findall("[-+]?[.]?[\d]+(?:,\d\d\d)*[\.]?\d*(?:[eE][-+]?\d+)?", fname)
    time,rank = float(rr[0]),int(rr[1])
    #print(fname,"->",time,rank)

    vtks_group[time].append(fname)


#Write into time-series XDMF format

vtks_time = {}
print('Time','Num_VTKFiles')
for time,vtks in vtks_group.items():
    vtks.sort(key=str.lower)
    print(time,len(vtks))

    blocks = pv.MultiBlock()
    for i in range(len(vtks)):
        blocks[f"Rank{i+1}"] = pv.read(vtks[i])
    vtk_union = blocks.combine()
    vtks_time[time] = vtk_union
    #print(time,vtk_union)
Times = list(vtks_time.keys())


data_t0=vtks_time[Times[0]]
data_meshio=pv2meshio(data_t0)

#Convert into a single time-series data with a shared mesh info
filename = 'solution_ngs.xdmf'
with meshio.xdmf.TimeSeriesWriter(filename) as writer:
    writer.write_points_cells(data_meshio.points, data_meshio.cells)
    for t in Times:
        p_data = vtks_time[t].point_arrays
        c_data = vtks_time[t].cell_arrays
        if(len(c_data)==0): c_data=None
        if(len(p_data)==0): p_data=None
        writer.write_data(t, point_data=p_data, cell_data=c_data)

print(f'Write time series mesh into {filename}')
	import netgen.meshing
	from netgen.geom2d import SplineGeometry


	from ngsolve import *

	comm = mpi_world
	rank = comm.rank
	np = comm.size

	do_vtk = True

	# viscosity
	nu = 0.001

	# timestepping parameters
	tau = 0.001
	tend = 3.0

	# mesh = Mesh("cylinder.vol")
	geo = SplineGeometry()
	geo.AddRectangle( (0, 0), (2, 0.41), bcs = ("wall", "outlet", "wall", "inlet"))
	geo.AddCircle ( (0.2, 0.2), r=0.05, leftdomain=0, rightdomain=1, bc="cyl")
	if rank==0:
	ngmesh = geo.GenerateMesh(maxh=0.08)
	ngmesh.Distribute(comm)
	else:
	ngmesh = netgen.meshing.Mesh.Receive(comm)
	ngmesh.SetGeometry(geo)

	mesh = Mesh(ngmesh)

	print(f"Rank{rank} NumEles={mesh.ne}")

	#does not work with mpi yet...
	#mesh.Curve(3)

	V = H1(mesh,order=3, dirichlet="wall\|cyl\|inlet")
	Q = H1(mesh,order=2)

	X = FESpace([V,V,Q])

	ux,uy,p = X.TrialFunction()
	vx,vy,q = X.TestFunction()

	div_u = grad(ux)[0]+grad(uy)[1]
	div_v = grad(vx)[0]+grad(vy)[1]

	stokes = nugrad(ux)grad(vx)+nugrad(uy)grad(vy)+div_uq+div_vp - 1e-10pq
	a = BilinearForm(X)
	a += SymbolicBFI(stokes)
	a.Assemble()

	# nothing here ...
	f = LinearForm(X)
	f.Assemble()

	# gridfunction for the solution
	gfu = GridFunction(X)

	# parabolic inflow at bc=1:
	uin = 1.54y(0.41-y)/(0.410.41)
	gfu.components[0].Set(uin, definedon=mesh.Boundaries("inlet"))

	velocity = CoefficientFunction(gfu.components[0:2])


	# solve Stokes problem for initial conditions:
	#inv_stokes = a.mat.Inverse(X.FreeDofs(), inverse="mumps")
	inv_stokes = a.mat.Inverse(X.FreeDofs(), inverse="masterinverse")
	res = f.vec.CreateVector()
	res.data = f.vec - a.mat*gfu.vec
	gfu.vec.data += inv_stokes * res

	# matrix for implicit Euler
	mstar = BilinearForm(X)
	mstar += SymbolicBFI(uxvx+uyvy + tau*stokes)
	mstar.Assemble()

	# inv = mstar.mat.Inverse(X.FreeDofs(), inverse="masterinverse")
	inv = mstar.mat.Inverse(X.FreeDofs(), inverse="masterinverse")

	# the non-linear term
	conv = BilinearForm(X, nonassemble = True)
	conv += SymbolicBFI( CoefficientFunction( (ux,uy) ) * (grad(ux)vx+grad(uy)vy) )

	t = 0
	vtk_interval = int(0.05/tau);

	import os
	#output_path = os.path.dirname(os.path.realpath(__file__)) + "/navierstokes_output"
	output_path = os.path.dirname(os.path.realpath(__file__)) + "/navierstokes_output"
	if rank==0 and not os.path.exists(output_path):
	os.mkdir(output_path)
	comm.Barrier() #wait until master has created the directory!!

	t=0.0

	if rank>0 and do_vtk:
	file_name = f"{output_path}/vtkout_time{t:.3f}_proc{rank} "
	print("rank "+str(rank)+", output to "+file_name)
	vtk = VTKOutput(ma=mesh,coefs=[velocity],names=["u"],filename=file_name,subdivision=1)
	vtk.Do()

	count = 1
	# implicit Euler/explicit Euler splitting method:
	with TaskManager():
	while t < tend:
	if rank==0:
	print ("t=", t)

	conv.Apply (gfu.vec, res)
	res.data += a.mat*gfu.vec
	gfu.vec.data -= tau * inv * res

	if rank>0 and count%vtk_interval==0 and do_vtk:
	file_name = f"{output_path}/vtkout_time{t:.3f}_proc{rank} "
	print("rank "+str(rank)+", output to "+file_name)
	vtk = VTKOutput(ma=mesh,coefs=[velocity],names=["u"],filename=file_name,subdivision=1)
	vtk.Do()
	count = count+1

	t = t + tau
	comm.Barrier()
	import numpy as np
	from numpy.core.shape_base import block
	import glob, os
	import re
	from collections import defaultdict

	import meshio
	import pyvista
	import pyvista as pv

	import argparse

	parser = argparse.ArgumentParser(description='This is a module to union parallel vtks for ngsolve.')
	parser.add_argument("--dir", default=".", help="Number of processor")
	args = parser.parse_args()
	print("Input Parameters=",args)


	try:
	from vtkmodules.vtkCommonCore import vtkVersion
	VTK9 = vtkVersion().GetVTKMajorVersion() >= 9
	except ImportError: # pragma: no cover
	VTK9 = False

	def pv2meshio(mesh):
	#save_meshio() in pyvista/utilities/fileio.py
	from meshio.vtk._vtk import vtk_to_meshio_type

	# Cast to pyvista.UnstructuredGrid
	if not isinstance(mesh, pv.UnstructuredGrid):
	mesh = mesh.cast_to_unstructured_grid()

	# Copy useful arrays to avoid repeated calls to properties
	vtk_offset = mesh.offset
	vtk_cells = mesh.cells
	vtk_cell_type = mesh.celltypes

	# Check that meshio supports all cell types in input mesh
	pixel_voxel = {8, 11} # Handle pixels and voxels
	for cell_type in np.unique(vtk_cell_type):
	if cell_type not in vtk_to_meshio_type.keys() and cell_type not in pixel_voxel:
	raise TypeError(f"meshio does not support VTK type {cell_type}.")

	# Get cells
	cells = []
	c = 0
	for offset, cell_type in zip(vtk_offset, vtk_cell_type):
	numnodes = vtk_cells[offset+c]
	if VTK9: # must offset by cell count
	cell = vtk_cells[offset+1+c:offset+1+c+numnodes]
	c += 1
	else:
	cell = vtk_cells[offset+1:offset+1+numnodes]
	cell = (
	cell if cell_type not in pixel_voxel
	else cell[[0, 1, 3, 2]] if cell_type == 8
	else cell[[0, 1, 3, 2, 4, 5, 7, 6]]
	)
	cell_type = cell_type if cell_type not in pixel_voxel else cell_type+1
	cell_type = (
	vtk_to_meshio_type[cell_type] if cell_type != 7
	else f"polygon{numnodes}"
	)

	if len(cells) > 0 and cells[-1][0] == cell_type:
	cells[-1][1].append(cell)
	else:
	cells.append((cell_type, [cell]))

	for k, c in enumerate(cells):
	cells[k] = (c[0], np.array(c[1]))


	return meshio.Mesh(np.array(mesh.points),cells)

	#Search vtk files in a given path
	os.chdir(args.dir)
	vtks=[f for f in glob.glob("*.vtk")]
	vtks.sort(key=str.lower)
	#print('VtkFiles=',vtks)

	#Parse filename and group vtk files based on time
	#ngsolve must using the following format
	#file_name = f"{output_path}/vtkout_time{t:.3f}_proc{rank} "

	vtks_group=defaultdict(list)
	for fname in vtks:
	#https://stackoverflow.com/questions/4289331/how-to-extract-numbers-from-a-string-in-python
	rr = re.findall("[-+]?[.]?[\d]+(?:,\d\d\d)[\.]?\d(?:[eE][-+]?\d+)?", fname)
	time,rank = float(rr[0]),int(rr[1])
	#print(fname,"->",time,rank)

	vtks_group[time].append(fname)


	#Write into time-series XDMF format

	vtks_time = {}
	print('Time','Num_VTKFiles')
	for time,vtks in vtks_group.items():
	vtks.sort(key=str.lower)
	print(time,len(vtks))

	blocks = pv.MultiBlock()
	for i in range(len(vtks)):
	blocks[f"Rank{i+1}"] = pv.read(vtks[i])
	vtk_union = blocks.combine()
	vtks_time[time] = vtk_union
	#print(time,vtk_union)
	Times = list(vtks_time.keys())


	data_t0=vtks_time[Times[0]]
	data_meshio=pv2meshio(data_t0)

	#Convert into a single time-series data with a shared mesh info
	filename = 'solution_ngs.xdmf'
	with meshio.xdmf.TimeSeriesWriter(filename) as writer:
	writer.write_points_cells(data_meshio.points, data_meshio.cells)
	for t in Times:
	p_data = vtks_time[t].point_arrays
	c_data = vtks_time[t].cell_arrays
	if(len(c_data)==0): c_data=None
	if(len(p_data)==0): p_data=None
	writer.write_data(t, point_data=p_data, cell_data=c_data)

	print(f'Write time series mesh into {filename}')