初始化仓库。

Author: zoziha

.gitignore

@@ -0,0 +1,4 @@
+build/
+*.o
+*.a
+*.x

Makefile

@@ -0,0 +1,25 @@
+# Fortran sph Makefile
+
+LIB = sph
+
+FC = gfortran
+FFLAGS = -g
+
+export LIB
+export FC
+export FFLAGS
+
+.PHONY: all clean test
+
+all:
+	$(MAKE) -f Makefile --directory=src
+	$(MAKE) -f Makefile --directory=app
+	$(MAKE) -f Makefile --directory=test
+
+test:
+	$(MAKE) -f Makefile --directory=test
+
+clean:
+	$(MAKE) -f Makefile clean --directory=src
+	$(MAKE) -f Makefile clean --directory=app
+	$(MAKE) -f Makefile clean --directory=test

README.md

@@ -0,0 +1,40 @@
+# 光滑粒子流体动力学（SPH）
+
+一份社区驱动的开源光滑粒子流体动力学（SPH）代码，起始代码版本源自课本《光滑粒子流体动力学--一种无网格粒子法》。
+
+| 项目 | 描述 |
+| --- | --- |
+| 版本 | 0.0.1 |
+| 许可证 | Public Domain |
+| 版权 | Copyright (c) 2021 SPH 贡献者 |
+
+## 开始
+
+### 获取代码
+
+```sh
+git clone https://github.com/zoziha/SPH.git
+cd SPH
+```
+
+### 使用Make构建代码
+
+本项目支持Make工具构建代码：
+
+```sh
+make
+```
+
+### 使用[fortran-lang/fpm](https://github.com/fortran-lang/fpm)构建代码
+
+FPM是社区驱动的Fortran语言的包管理器和代码构建工具，适用于c/c++/fortran代码的构建。  
+你可以通过提供的`fpm.toml`构建代码：
+
+```sh
+fpm build
+fpm run
+```
+
+## 链接
+
++ [spheric/SPH Codes](https://spheric-sph.org/sph-projects-and-codes)

app/Makefile

@@ -0,0 +1,6 @@
+all: sph
+
+# Orginal test
+sph: sph.o 
+	$(FC) $(FFLAGS) $< -L../src -l$(LIB) -I../src -o $@.x
+	cd .. && time ./app/sph.x

app/sph.f

@@ -0,0 +1,52 @@
+      program SPH
+
+c----------------------------------------------------------------------
+c     This is a three dimensional SPH code. the followings are the 
+c     basic parameters needed in this code or calculated by this code
+
+c     mass-- mass of particles                                      [in]
+c     ntotal-- total particle number ues                            [in]
+c     dt--- Time step used in the time integration                  [in]
+c     itype-- types of particles                                    [in]
+c     x-- coordinates of particles                              [in/out]
+c     vx-- velocities of particles                              [in/out]
+c     rho-- densities of particles                              [in/out]
+c     p-- pressure  of particles                                [in/out]
+c     u-- internal energy of particles                          [in/out]
+c     hsml-- smoothing lengths of particles                     [in/out]
+c     c-- sound velocity of particles                              [out]
+c     s-- entropy of particles                                     [out]
+c     e-- total energy of particles                                [out]
+
+      implicit none     
+      include '../src/param.inc'
+
+      integer ntotal, itype(maxn), maxtimestep, d, m, i, yesorno      
+      double precision x(dim,maxn), vx(dim, maxn), mass(maxn),rho(maxn),
+     &     p(maxn), u(maxn), c(maxn), s(maxn), e(maxn), hsml(maxn), dt
+      double precision s1, s2
+
+      call time_print
+      call time_elapsed(s1)      
+
+      if (shocktube)   dt = 0.005
+      if (shearcavity) dt = 5.e-5
+      call input(x, vx, mass, rho, p, u, itype, hsml, ntotal)     
+ 1    write(*,*)'  ***************************************************' 
+      write(*,*)'          Please input the maximal time steps '
+      write(*,*)'  ***************************************************'
+      read(*,*) maxtimestep      
+      call time_integration(x, vx, mass, rho, p, u, c, s, e, itype, 
+     &     hsml, ntotal, maxtimestep, dt )
+      call output(x, vx, mass, rho, p, u, c, itype, hsml, ntotal)      
+      write(*,*)'  ***************************************************'
+      write(*,*) 'Are you going to run more time steps ? (0=No, 1=yes)'
+      write(*,*)'  ***************************************************'
+      read (*,*) yesorno     
+      if (yesorno.ne.0) go to 1
+      call time_print
+      call time_elapsed(s2)      
+      write (*,*)'        Elapsed CPU time = ', s2-s1
+	  write (*,*)'all finish!'
+                           
+      end

fpm.toml

@@ -0,0 +1,14 @@
+name = "sph_code"
+
+[build]
+auto-executables = false
+auto-tests = false
+auto-examples = false
+
+[library]
+include-dir = "src"
+
+[[executable]]
+name = "sph"
+source-dir = "app"
+main = "sph.f"

src/Makefile

@@ -0,0 +1,41 @@
+.SUFFIXES: .f90 .o
+
+SRCF = \
+	art_heat.f \
+	art_visc.f \
+	av_vel.f \
+	density.f \
+	direct_find.f \
+	eos.f \
+	external_force.f \
+	grid_geom.f \
+	hsml.f \
+	init_grid.f \
+	input.f \
+	internal_force.f \
+	kernel.f \
+	link_list.f \
+	output.f \
+	single_step.f \
+	time_integration.f \
+	virt_part.f \
+	viscosity.f
+
+SRCF90 = \
+	time_elapsed.f90 \
+	time_print.f90
+
+OBJF   := $(SRCF:.f=.o)
+OBJF90 := $(SRCF90:.f90=.o)
+
+lib$(LIB).a: $(OBJF) $(OBJF90)
+	ar -rcs lib$(LIB).a $(OBJF) $(OBJF90)
+
+clean:
+	rm -f -r *.o *.a *.so *.mod *.smod
+
+%.o: %.f90
+	$(FC) $(FFLAGS) -c $<
+
+time_elapsed.o: time_elapsed.f90
+time_print.o  : time_print.f90

src/art_heat.f

@@ -0,0 +1,83 @@
+      subroutine art_heat(ntotal,hsml,mass,x,vx,niac,rho,u,
+     &           c,pair_i,pair_j,w,dwdx,dedt)
+
+c----------------------------------------------------------------------
+c     Subroutine to calculate the artificial heat(Fulk, 1994, p, a-17) 
+c     See Equ.(4.74)  
+
+c     ntotal : Number of particles                                  [in]
+c     hsml   : Smoothing Length                                     [in]
+c     mass   : Particle masses                                      [in]
+c     x      : Coordinates of all particles                         [in]
+c     vx     : Velocities of all particles                          [in]
+c     rho    : Density                                              [in]
+c     u      : specific internal energy                             [in]
+c     c      : Sound veolcity                                       [in]
+c     niac   : Number of interaction pairs                          [in]
+c     pair_i : List of first partner of interaction pair            [in]
+c     pair_j : List of second partner of interaction pair           [in]
+c     w      : Kernel for all interaction pairs                     [in]
+c     dwdx   : Derivative of kernel with respect to x, y and z      [in]
+c     dedt   : produced artificial heat, adding to energy Eq.      [out]
+
+      implicit none
+      include 'param.inc'
+      
+      integer ntotal,niac,pair_i(max_interaction),
+     &        pair_j(max_interaction)
+      double precision hsml(maxn), mass(maxn), x(dim,maxn),vx(dim,maxn),
+     &                 rho(maxn), u(maxn), c(maxn),w(max_interaction),
+     &                 dwdx(dim,max_interaction), dedt(maxn)
+      integer i,j,k,d
+      double precision dx, dvx(dim), vr, rr, h, mc, mrho, mhsml, 
+     &                 vcc(maxn), hvcc, mui, muj, muij, rdwdx, g1,g2
+      
+
+c---  Parameter for the artificial heat conduction:
+     
+      g1=0.1
+      g2=1.0      
+      do i=1,ntotal
+        vcc(i) = 0.e0         
+        dedt(i) = 0.e0
+      enddo
+    
+      do k=1,niac
+        i = pair_i(k)
+        j = pair_j(k)
+        do d=1,dim
+          dvx(d) = vx(d,j) - vx(d,i) 
+        enddo        
+        hvcc = dvx(1)*dwdx(1,k)
+        do d=2,dim
+          hvcc = hvcc + dvx(d)*dwdx(d,k)
+        enddo    
+         vcc(i) = vcc(i) + mass(j)*hvcc/rho(j)
+         vcc(j) = vcc(j) + mass(i)*hvcc/rho(i)
+      enddo  
+   
+      do k=1,niac
+        i = pair_i(k)
+        j = pair_j(k)
+        mhsml= (hsml(i)+hsml(j))/2.
+        mrho = 0.5e0*(rho(i) + rho(j))                          
+        rr = 0.e0
+        rdwdx = 0.e0
+        do d=1,dim
+          dx = x(d,i) -  x(d,j)
+          rr = rr + dx*dx
+          rdwdx  = rdwdx + dx*dwdx(d,k)            
+        enddo             
+        mui=g1*hsml(i)*c(i) + g2*hsml(i)**2*(abs(vcc(i))-vcc(i))
+        muj=g1*hsml(j)*c(j) + g2*hsml(j)**2*(abs(vcc(j))-vcc(j))
+        muij= 0.5*(mui+muj)     
+        h = muij/(mrho*(rr+0.01*mhsml**2))*rdwdx
+        dedt(i) = dedt(i) + mass(j)*h*(u(i)-u(j))
+        dedt(j) = dedt(j) + mass(i)*h*(u(j)-u(i))
+      enddo
+
+      do i=1,ntotal
+        dedt(i) = 2.0e0*dedt(i)          
+      enddo
+
+      end

src/art_visc.f

@@ -0,0 +1,99 @@
+      subroutine art_visc(ntotal,hsml,mass,x,vx,niac,rho,c,
+     &           pair_i,pair_j,w,dwdx,dvxdt,dedt)
+
+c----------------------------------------------------------------------
+c     Subroutine to calculate the artificial viscosity (Monaghan, 1992) 
+c     See Equ.(4.66) Equ.(4.62)
+
+c     ntotal : Number of particles (including virtual particles)    [in]
+c     hsml   : Smoothing Length                                     [in]
+c     mass   : Particle masses                                      [in]
+c     x      : Coordinates of all particles                         [in]
+c     vx     : Velocities of all particles                          [in]
+c     niac   : Number of interaction pairs                          [in]
+c     rho    : Density                                              [in]
+c     c      : Temperature                                          [in]
+c     pair_i : List of first partner of interaction pair            [in]
+c     pair_j : List of second partner of interaction pair           [in]
+c     w      : Kernel for all interaction pairs                     [in]
+c     dwdx   : Derivative of kernel with respect to x, y and z      [in]
+c     dvxdt  : Acceleration with respect to x, y and z             [out] 
+c     dedt   : Change of specific internal energy                  [out]
+ 
+      implicit none
+      include 'param.inc'
+      
+      integer ntotal, niac, pair_i(max_interaction),
+     &        pair_j(max_interaction)            
+      double precision hsml(maxn), mass(maxn), x(dim,maxn),vx(dim,maxn),
+     &       rho(maxn), c(maxn), w(max_interaction),
+     &       dwdx(dim,max_interaction), dvxdt(dim,maxn), dedt(maxn)
+      integer i,j,k,d
+      double precision dx, dvx(dim), alpha, beta, etq, piv,
+     &       muv, vr, rr, h, mc, mrho, mhsml
+     
+c     Parameter for the artificial viscosity:
+c     Shear viscosity
+      parameter( alpha = 1.e0   )
+     
+c     Bulk viscosity
+      parameter( beta  = 1.e0  ) 
+      
+c     Parameter to avoid singularities
+      parameter( etq   = 0.1e0 )
+           
+      do i=1,ntotal
+        do d=1,dim
+          dvxdt(d,i) = 0.e0
+        enddo
+        dedt(i) = 0.e0
+      enddo   
+     
+c     Calculate SPH sum for artificial viscosity
+      
+      do k=1,niac
+        i = pair_i(k)
+        j = pair_j(k)
+        mhsml= (hsml(i)+hsml(j))/2.
+        vr = 0.e0
+        rr = 0.e0
+        do d=1,dim
+          dvx(d) = vx(d,i) - vx(d,j)
+          dx     =  x(d,i) -  x(d,j)
+          vr     = vr + dvx(d)*dx
+          rr     = rr + dx*dx
+        enddo
+
+c     Artificial viscous force only if v_ij * r_ij < 0
+
+        if (vr.lt.0.e0) then
+
+c     Calculate muv_ij = hsml v_ij * r_ij / ( r_ij^2 + hsml^2 etq^2 )
+            
+          muv = mhsml*vr/(rr + mhsml*mhsml*etq*etq)
+          
+c     Calculate PIv_ij = (-alpha muv_ij c_ij + beta muv_ij^2) / rho_ij
+
+          mc   = 0.5e0*(c(i) + c(j))
+          mrho = 0.5e0*(rho(i) + rho(j))
+          piv  = (beta*muv - alpha*mc)*muv/mrho              
+
+c     Calculate SPH sum for artificial viscous force
+
+          do d=1,dim
+            h = -piv*dwdx(d,k)
+            dvxdt(d,i) = dvxdt(d,i) + mass(j)*h
+            dvxdt(d,j) = dvxdt(d,j) - mass(i)*h
+            dedt(i) = dedt(i) - mass(j)*dvx(d)*h
+            dedt(j) = dedt(j) - mass(i)*dvx(d)*h
+          enddo
+        endif
+      enddo
+
+c     Change of specific internal energy:
+
+      do i=1,ntotal
+         dedt(i) = 0.5e0*dedt(i)        
+      enddo
+
+      end