From cb9b398a0d4ff55b2a646c81c3cebd2cfd7c36c8 Mon Sep 17 00:00:00 2001 From: Manas Rachh Date: Wed, 21 Feb 2024 12:15:42 -0500 Subject: [PATCH] adding stokes fmm python interfaces --- docs/conf.py | 1 - python/fmm2dpy/__init__.py | 2 +- python/fmm2dpy/fmm2d.py | 337 +++++++++++++++++++++++++++++++++++++ python/setup.py | 14 +- python/stfmm2d_example.py | 36 ++++ src/stokes/stfmm2d.f | 23 ++- 6 files changed, 402 insertions(+), 11 deletions(-) create mode 100644 python/stfmm2d_example.py diff --git a/docs/conf.py b/docs/conf.py index 8d3f29d..c6128d5 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -37,7 +37,6 @@ 'sphinx.ext.todo', 'sphinx.ext.napoleon', 'sphinx.ext.mathjax', -# 'sphinx.ext.autosectionlabel', # needs v 1.4; can :ref: other files w/o this; removed 7/29/18 'texext', # 'sphinxcontrib.bibtex', ] diff --git a/python/fmm2dpy/__init__.py b/python/fmm2dpy/__init__.py index f97a7d1..4ccde61 100644 --- a/python/fmm2dpy/__init__.py +++ b/python/fmm2dpy/__init__.py @@ -1 +1 @@ -from .fmm2d import hfmm2d,rfmm2d,lfmm2d,cfmm2d,bhfmm2d,h2ddir,r2ddir,l2ddir,c2ddir,bh2ddir,comperr,Output +from .fmm2d import hfmm2d,rfmm2d,lfmm2d,cfmm2d,bhfmm2d,stfmm2d,h2ddir,r2ddir,l2ddir,c2ddir,bh2ddir,st2ddir,comperr,Output diff --git a/python/fmm2dpy/fmm2d.py b/python/fmm2dpy/fmm2d.py index 5d9b977..a6f4a76 100644 --- a/python/fmm2dpy/fmm2d.py +++ b/python/fmm2dpy/fmm2d.py @@ -1,6 +1,7 @@ from . import hfmm2d_fortran as hfmm from . import lfmm2d_fortran as lfmm from . import bhfmm2d_fortran as bhfmm +from . import stfmm2d_fortran as stfmm import numpy as np import numpy.linalg as la @@ -12,6 +13,8 @@ class Output(): pottarg = None gradtarg = None hesstarg = None + pre = None + pretarg = None ier = 0 def hfmm2d(*,eps,zk,sources,charges=None,dipstr=None,dipvec=None, @@ -1094,6 +1097,169 @@ def bhfmm2d(*,eps,sources,charges=None,dipoles=None, return out +def stfmm2d(*,eps,sources,stoklet=None,strslet=None,strsvec=None,targets=None,ifppreg=0,ifppregtarg=0,nd=1): + r""" + Stokes FMM in R^{2}: evaluate all pairwise particle + interactions (ignoring self-interactions) and + interactions with targs. + + This routine computes sums of the form + + u(x) = sum_m G_{ij}(x,y^{(m)}) sigma^{(m)}_j + + sum_m T_{ijk}(x,y^{(m)}) mu^{(m)}_j nu^{(m)}_k + + where sigma^{(m)} is the Stokeslet charge, mu^{(m)} is the + stresslet charge, and nu^{(m)} is the stresslet orientation + (note that each of these is a 2 vector per source point y^{(m)}). + For x a source point, the self-interaction in the sum is omitted. + + Optionally, the associated pressure p(x) and gradient grad u(x) + are returned + + p(x) = sum_m P_j(x,y^m) sigma^{(m)}_j + + sum_m T_{ijk}(x,y^{(m)}) PI_{jk} mu^{(m)}_j nu^{(m)}_k + + grad u(x) = grad[sum_m G_{ij}(x,y^m) sigma^{(m)}_j + + sum_m T_{ijk}(x,y^{(m)}) mu^{(m)}_j nu^{(m)}_k] + + where G_{ij} is the Stokelet, P_{j} its associated pressure + tensor, T_{ijk} is the type I stresslet, and PI_{jk} is its + associated pressure tensor. These kernels are given by + + G_{ij} = (r_{i} r_{j})/(2 r^2) - \delta_{ij} \log(r)/2 + P_{j} = r_{j}/r^2 + T_{ijk} = - 2 r_{i} r_{j} r_{k}/r^4 + PI_{jk} = -\delta_{jk}/r^2 + 2r_{j} r_{k}/r^{4} + + where for a source y, and target x, r_{i} = x_{i}-y_{i}, + and r = \sqrt{r_{1}^2 + r_{2}^2} + + Args: + eps: float + precision requested + sources: float(2,n) + source locations + stoklet: float(nd,2,n) or float(2,n) + Stokeslet charge strengths (sigma vectors above) + strslet: float(nd,2,n) or float(2,n) + stresslet strengths (mu vectors above) + strsvec: float(nd,2,n) or float(2,n) + stresslet orientations (nu vectors above) + targets: float(2,nt) + target locations (x) + + ifppreg: integer + flag for evaluating potential, gradient, and pressure + at the sources + ifppreg = 1, only potential + ifppreg = 2, potential and pressure + ifppreg = 3, potential, pressure, and gradient + + ifppregtarg: integer + flag for evaluating potential, gradient, and pressure + at the targets + ifppregtarg = 1, only potential + ifppregtarg = 2, potential and pressure + ifppregtarg = 3, potential, pressure, and gradient + + nd: integer + number of densities + + Returns: + out.pot: velocity at source locations if requested + out.pre: pressure at source locations if requested + out.grad: gradient of velocity at source locations if requested + out.pottarg: velocity at target locations if requested + out.pretarg: pressure at target locations if requested + out.gradtarg: gradient of velocity at target locations if requested + + Example: + see stfmmexample.py + + r""" + out = Output() + + if(ifppreg == 0 and ifppregtarg == 0): + print("Nothing to compute, set either ifppreg or ifppregtarg to non-zero") + return out + + if(stoklet is None and strslet is None and strsvec is None): + print("Nothing to compute, set either stoklet or strslet+strsvec to non-None") + return out + + if(strslet is not None and strsvec is None): + print("strslet and strsvec mush be both None or both not None") + return out + if(strslet is None and strsvec is not None): + print("strslet and strsvec mush be both None or both not None") + return out + + assert sources.shape[0] == 2, "The first dimension of sources must be 2" + if(np.size(np.shape(sources))==2): + ns = sources.shape[1] + if(np.size(np.shape(sources))==1): + ns = 1 + if(targets is not None): + assert targets.shape[0] == 2, "The first dimension of targets must be 2" + if(np.size(np.shape(targets))==2): + nt = targets.shape[1] + if(np.size(np.shape(targets))==1): + nt = 1 + else: + targets = np.zeros([2,0],dtype='double') + nt = 0 + + ifstoklet = 0 + ifstrslet = 0 + + if(stoklet is not None): + if(nd == 1 and ns>1): + assert stoklet.shape[0] == 2 and stoklet.shape[1] == ns, "stoklet vectors must be of shape [2,number of sources]" + if(nd == 1 and ns==1): + assert stoklet.shape[0] == 2, "stoklet vectors must be of shape [2,number of sources]" + if(nd>1): + assert stoklet.shape[0] == nd and stoklet.shape[1] == 2 and stoklet.shape[2] == ns, "stoklet vectors must be of shape [nd,2,ns] where nd is number of densities, and ns is number of sources" + stoklet = stoklet.reshape([nd,2,ns]) + ifstoklet = 1 + else: + stoklet = np.zeros([nd,2,ns],dtype='double') + + if(strslet is not None and strsvec is not None): + if(nd == 1 and ns>1): + assert strslet.shape[0] == 2 and strslet.shape[1] == ns, "strslet vectors must be of shape [2,number of sources]" + assert strsvec.shape[0] == 2 and strsvec.shape[1] == ns, "strsvec vectors must be of shape [2,number of sources]" + if(nd == 1 and ns==1): + assert strslet.shape[0] == 2, "strslet vectors must be of shape [2,number of sources]" + assert strsvec.shape[0] == 2, "strsvec vectors must be of shape [2,number of sources]" + if(nd>1): + assert strslet.shape[0] == nd and strslet.shape[1] == 2 and strslet.shape[2] == ns, "strslet vectors must be of shape [nd,2,ns] where nd is number of densities, and ns is number of sources" + assert strsvec.shape[0] == nd and strsvec.shape[1] == 2 and strsvec.shape[2] == ns, "strsvec vectors must be of shape [nd,2,ns] where nd is number of densities, and ns is number of sources" + strslet = strslet.reshape([nd,2,ns]) + strsvec = strsvec.reshape([nd,2,ns]) + ifstrslet = 1 + else: + strslet = np.zeros([nd,2,ns],dtype='double') + strsvec = np.zeros([nd,2,ns],dtype='double') + + out.pot,out.pre,out.grad,out.pottarg,out.pretarg,out.gradtarg,out.ier = stfmm.stfmm2d(eps,sources,ifstoklet,stoklet,ifstrslet,strslet,strsvec,ifppreg,targets,ifppregtarg,nd,ns,nt) + + if(ifppreg < 3): + out.grad = None + if(ifppregtarg < 3): + out.gradtarg = None + + if(ifppreg < 2): + out.pre = None + if(ifppregtarg < 2): + out.pretarg = None + + if(ifppreg < 1): + out.pot = None + if(ifppregtarg < 1): + out.pottarg = None + + return out + def h2ddir(*,zk,sources,targets,charges=None,dipstr=None,dipvec=None, pgt=0,nd=1,thresh=1e-16): r""" @@ -1689,6 +1855,177 @@ def bh2ddir(*,sources,targets,charges=None,dipoles=None, return out +def st2ddir(*,eps,sources,stoklet=None,strslet=None,strsvec=None,targets=None,ifppreg=0,ifppregtarg=0,nd=1,thresh=1e-16): + r""" + This subroutine evaluates all pairwise particle + interactions (ignoring self-interactions) and + interactions with targs. + + This routine computes sums of the form + + u(x) = sum_m G_{ij}(x,y^{(m)}) sigma^{(m)}_j + + sum_m T_{ijk}(x,y^{(m)}) mu^{(m)}_j nu^{(m)}_k + + where sigma^{(m)} is the Stokeslet charge, mu^{(m)} is the + stresslet charge, and nu^{(m)} is the stresslet orientation + (note that each of these is a 2 vector per source point y^{(m)}). + For x a source point, the self-interaction in the sum is omitted. + + Optionally, the associated pressure p(x) and gradient grad u(x) + are returned + + p(x) = sum_m P_j(x,y^m) sigma^{(m)}_j + + sum_m T_{ijk}(x,y^{(m)}) PI_{jk} mu^{(m)}_j nu^{(m)}_k + + grad u(x) = grad[sum_m G_{ij}(x,y^m) sigma^{(m)}_j + + sum_m T_{ijk}(x,y^{(m)}) mu^{(m)}_j nu^{(m)}_k] + + where G_{ij} is the Stokelet, P_{j} its associated pressure + tensor, T_{ijk} is the type I stresslet, and PI_{jk} is its + associated pressure tensor. These kernels are given by + + G_{ij} = (r_{i} r_{j})/(2 r^2) - \delta_{ij} \log(r)/2 + P_{j} = r_{j}/r^2 + T_{ijk} = - 2 r_{i} r_{j} r_{k}/r^4 + PI_{jk} = -\delta_{jk}/r^2 + 2r_{j} r_{k}/r^{4} + + where for a source y, and target x, r_{i} = x_{i}-y_{i}, + and r = \sqrt{r_{1}^2 + r_{2}^2} + + Args: + eps: float + precision requested + sources: float(2,n) + source locations + stoklet: float(nd,2,n) or float(2,n) + Stokeslet charge strengths (sigma vectors above) + strslet: float(nd,2,n) or float(2,n) + stresslet strengths (mu vectors above) + strsvec: float(nd,2,n) or float(2,n) + stresslet orientations (nu vectors above) + targets: float(2,nt) + target locations (x) + + ifppreg: integer + flag for evaluating potential, gradient, and pressure + at the sources + ifppreg = 1, only potential + ifppreg = 2, potential and pressure + ifppreg = 3, potential, pressure, and gradient + + ifppregtarg: integer + flag for evaluating potential, gradient, and pressure + at the targets + ifppregtarg = 1, only potential + ifppregtarg = 2, potential and pressure + ifppregtarg = 3, potential, pressure, and gradient + + nd: integer + number of densities + + thresh: contribution of source x_i, at location x ignored if |x-x_i|<=thresh + + Returns: + out.pot: velocity at source locations if requested + out.pre: pressure at source locations if requested + out.grad: gradient of velocity at source locations if requested + out.pottarg: velocity at target locations if requested + out.pretarg: pressure at target locations if requested + out.gradtarg: gradient of velocity at target locations if requested + + Example: + see stfmmexample.py + + r""" + out = Output() + + if(ifppreg == 0 and ifppregtarg == 0): + print("Nothing to compute, set either ifppreg or ifppregtarg to non-zero") + return out + + if(stoklet == None and strslet == None and strsvec == None): + print("Nothing to compute, set either stoklet or strslet+strsvec to non-None") + return out + + if(strslet is not None and strsvec is None): + print("strslet and strsvec mush be both None or both not None") + return out + if(strslet is None and strsvec is not None): + print("strslet and strsvec mush be both None or both not None") + return out + + assert sources.shape[0] == 2, "The first dimension of sources must be 2" + if(np.size(np.shape(sources))==2): + ns = sources.shape[1] + if(np.size(np.shape(sources))==1): + ns = 1 + if(targets is not None): + assert targets.shape[0] == 2, "The first dimension of targets must be 2" + if(np.size(np.shape(targets))==2): + nt = targets.shape[1] + if(np.size(np.shape(targets))==1): + nt = 1 + else: + targets = np.zeros([2,0],dtype='double') + nt = 0 + + ifstoklet = 0 + ifstrslet = 0 + + if(stoklet is not None): + if(nd == 1 and ns>1): + assert stoklet.shape[0] == 2 and stoklet.shape[1] == ns, "stoklet vectors must be of shape [2,number of sources]" + if(nd == 1 and ns==1): + assert stoklet.shape[0] == 2, "stoklet vectors must be of shape [2,number of sources]" + if(nd>1): + assert stoklet.shape[0] == nd and stoklet.shape[1] == 2 and stoklet.shape[2] == ns, "stoklet vectors must be of shape [nd,2,ns] where nd is number of densities, and ns is number of sources" + stoklet = stoklet.reshape([nd,2,ns]) + ifstoklet = 1 + else: + stoklet = np.zeros([nd,2,ns],dtype='double') + + if(strslet is not None and strsvec is not None): + if(nd == 1 and ns>1): + assert strslet.shape[0] == 2 and strslet.shape[1] == ns, "strslet vectors must be of shape [2,number of sources]" + assert strsvec.shape[0] == 2 and strsvec.shape[1] == ns, "strsvec vectors must be of shape [2,number of sources]" + if(nd == 1 and ns==1): + assert strslet.shape[0] == 2, "strslet vectors must be of shape [2,number of sources]" + assert strsvec.shape[0] == 2, "strsvec vectors must be of shape [2,number of sources]" + if(nd>1): + assert strslet.shape[0] == nd and strslet.shape[1] == 2 and strslet.shape[2] == ns, "strslet vectors must be of shape [nd,2,ns] where nd is number of densities, and ns is number of sources" + assert strsvec.shape[0] == nd and strsvec.shape[1] == 2 and strsvec.shape[2] == ns, "strsvec vectors must be of shape [nd,2,ns] where nd is number of densities, and ns is number of sources" + strslet = strslet.reshape([nd,2,ns]) + strsvec = strsvec.reshape([nd,2,ns]) + ifstrslet = 1 + else: + strslet = np.zeros([nd,2,ns],dtype='double') + strsvec = np.zeros([nd,2,ns],dtype='double') + + if(ifstoklet == 1 and ifstrslet == 0): + out.pot,out.pre,out.grad = stfmm.st2ddirectstokg(sources,stoklet,sources,thresh,nd,ns,nt) + out.pottarg,out.pretarg,out.gradtarg = stfmm.st2ddirectstokg(sources,stoklet,targets,thresh,nd,ns,nt) + else: + out.pot,out.pre,out.grad = stfmm.st2ddirectstokstrsg(sources,stoklet,1,strslet,strsvec,sources,thresh,nd,ns,nt) + out.pottarg,out.pretarg,out.gradtarg = stfmm.st2ddirectstokstrsg(sources,stoklet,1,strslet,strsvec,targets,thresh,nd,ns,nt) + + if(ifppreg < 3): + out.grad = None + if(ifppregtarg < 3): + out.gradtarg = None + + if(ifppreg < 2): + out.pre = None + if(ifppregtarg < 2): + out.pretarg = None + + if(ifppreg < 1): + out.pot = None + if(ifppregtarg < 1): + out.pottarg = None + + return out + + def comperr(*,ntest,out,outex,pg=0,pgt=0,nd=1,cauchy=0,bh=0): r = 0 diff --git a/python/setup.py b/python/setup.py index 874ed9f..3211b8c 100644 --- a/python/setup.py +++ b/python/setup.py @@ -13,6 +13,7 @@ list_helm=['hfmm2dwrap.f','hfmm2dwrap_vec.f','helmkernels2d.f'] list_lap=['rfmm2dwrap.f','rfmm2dwrap_vec.f','rlapkernels2d.f','lfmm2dwrap.f','lfmm2dwrap_vec.f','lapkernels2d.f','cfmm2dwrap.f','cfmm2dwrap_vec.f','cauchykernels2d.f'] list_bh=['bhfmm2dwrap.f','bhkernels2d.f'] +list_stk=['stfmm2d.f','stokkernels2d.f'] list_common=[] FLIBS = os.getenv('FMM_FLIBS') @@ -94,24 +95,31 @@ extra_link_args=FLIBS ) +ext_st = Extension( + name='fmm2dpy.stfmm2d_fortran', + sources=['../src/stokes/stfmm2d.f']+['../src/stokes/stokkernels2d.f'], + f2py_options=['only:']+['stfmm2d']+['st2ddirectstokg']+['st2ddirectstokstrsg']+[':'], + extra_link_args=FLIBS +) + ## TODO: fill in the info below setup( name=pkg_name, python_requires='>=3.0.0', - version="0.0.5", + version="1.1.0", author="Travis Askham, Zydrunas Gimbutas, Leslie Greengard, Libin Lu, Michael O'Neil, Manas Rachh, and Vladimir Rokhlin", author_email="mrachh@flatironinstitute.org", description="This pacakge contains basic routines for Laplace and Helmholtz fast multipole methods in two dimensions", long_description=open('../README.md').read(), long_description_content_type='text/markdown', - url="", + url="https://github.com/flatironinstitute/fmm2d", packages=['fmm2dpy'], install_requires=[ "numpy", "pytest" ], - ext_modules=[ext_helm,ext_lap,ext_bh], + ext_modules=[ext_helm,ext_lap,ext_bh,ext_st], classifiers=[ "Programming Language :: Python :: 3", "License :: OSI Approved :: Apache Software License", diff --git a/python/stfmm2d_example.py b/python/stfmm2d_example.py new file mode 100644 index 0000000..1c396e8 --- /dev/null +++ b/python/stfmm2d_example.py @@ -0,0 +1,36 @@ +#!/usr/bin/env python + +import fmm2dpy as fmm +import numpy as np + + +# +# This is a sample code to demonstrate how to use +# the fmm libraries +# + +# sample with one density, sources to sources, +# stokeslet interactions, and potential+pressure+gradient +# +n = 20000 +nd = 1 +sources = np.random.uniform(0,1,(2,n)) +eps = 10**(-5) + +stoklet = np.random.uniform(0,1,(2,n)) +out = fmm.stfmm2d(eps=eps,sources=sources,stoklet=stoklet,ifppreg=3) + + +# sample with a vector of densities, sources to +# sources and targets, stokeslet + stresslet interactions, +# potential + pressure + +nd = 2 +nt = 1870 +targ = np.random.uniform(0,1,(2,nt)) +stoklet = np.random.uniform(0,1,(nd,2,n)) +strsvec = np.random.uniform(0,1,(nd,2,n)) +strslet = np.random.uniform(0,1,(nd,2,n)) +out = fmm.stfmm2d(eps=eps,sources=sources,stoklet=stoklet,\ + strslet=strslet,strsvec=strsvec,\ + targets=targ,nd=nd,ifppreg=2,ifppregtarg=2) diff --git a/src/stokes/stfmm2d.f b/src/stokes/stfmm2d.f index f8f9d3a..d5fda62 100644 --- a/src/stokes/stfmm2d.f +++ b/src/stokes/stfmm2d.f @@ -29,6 +29,16 @@ subroutine stfmm2d(nd, eps, $ ifstoklet, stoklet, ifstrslet, strslet, strsvec, $ ifppreg, pot, pre, grad, ntarg, targ, $ ifppregtarg, pottarg, pretarg, gradtarg,ier) + +cf2py intent(in) nd,eps +cf2py intent(in) nsource,source +cf2py intent(in) ifstoklet,stoklet +cf2py intent(in) ifstrslet,strslet,strsvec +cf2py intent(in) ifppreg,ifppregtarg +cf2py intent(in) ntarg,targ +cf2py intent(out) pot,pre,grad +cf2py intent(out) pottarg,pretarg,gradtarg +cf2py intent(out) ier c c Stokes FMM in R^{3}: evaluate all pairwise particle c interactions (ignoring self-interactions) and @@ -154,12 +164,13 @@ subroutine stfmm2d(nd, eps, integer nd, ifstoklet, ifstrslet, ntarg double precision eps integer nsource, ifppreg, ifppregtarg - double precision source(2, *), targ(2, *) - double precision stoklet(nd, 2, *), strslet(nd, 2, *) - double precision strsvec(nd, 2, *) - double precision pot(nd, 2, *), pre(nd,*), grad(nd, 2, 2, *) - double precision pottarg(nd, 2, *), pretarg(nd,*), - 1 gradtarg(nd, 2, 2, *) + double precision source(2, nsource), targ(2, ntarg) + double precision stoklet(nd, 2, nsource), strslet(nd, 2, nsource) + double precision strsvec(nd, 2, nsource) + double precision pot(nd, 2, nsource), pre(nd, nsource) + double precision grad(nd, 2, 2, nsource) + double precision pottarg(nd, 2, ntarg), pretarg(nd,ntarg), + 1 gradtarg(nd, 2, 2, ntarg) c c CONTINUE FROM HERE c