|
|
実行時内部エラーが発生しました。
access violation エラーです。
decimalBASIC 6.6.3.1 / MacOS 10.7.5 / MACmini (intel core-i7 2.7GHz 4GB)
再起動後に実行してもエラーが再現します。
下記のプログラムです。
コメントアウトして発生場所を絞り込むと、
SUB setInitialCondition で発生しているものと推定できます。
!set density
FOR i=0 TO NNp-1
LET s = 0.0
FOR j=0 TO NNp-1
LET s = s + mass(j)*kernel(distance(i,j), (hh(i)+hh(j))/2.0)
NEXT j
LET density(i) = s
NEXT i
の部分を下記のように変更するとエラーは回避できそうですが
!set density
FOR i=0 TO NNp-1
LET s = 0.0
FOR j=0 TO NNp-1
LET r = distance(i,j)
LET s = s + mass(j)*kernel(r, (hh(i)+hh(j))/2.0)
NEXT j
LET density(i) = s
NEXT i
まだ未完ですが、全体のプログラムを載せます。
-------------
!
! ========= smoothed particle hydrodynamics 2D ==========
!
! 0xxgasSPH2D.bas
!
! Copyright(C) 2017 Mitsuru Ikeuchi
!
! ver 0.0.0 2017.06.22 created
!
OPTION ARITHMETIC NATIVE
DECLARE EXTERNAL SUB sph2d.setInitialCondition, sph2d.timeEvolution, sph2d.drawParticles
CALL setInitialCondition
LET t0 = TIME
FOR it=1 TO 1000
!CALL timeEvolution(10)
!CALL drawParticles
NEXT it
PRINT TIME - t0
END
! ---------- smoothed particle hydrodynamics 2D ----------
!
! - particle base Lagrangian method
! Monaghan; "Smoothed Particle Hydrodynamics"
! Annu. Rev. Astron. Astrophys.1992. 30:543-74
!
! W(x-xi,h) : kernel weight function (q=|x-xi|/h)
! = aKernel*(1-1.5*q^2+0.75*q^3) (q<1)
! = aKernel*0.25*(2-q)^3 (1<=q<2)
! = 0 (q>=2)
!
! aKernel = 2/3/h (1D)
! = 10/(7Pi)/h^2 (2D)
! = 1/Pi/h^3 (3D)
!
! f(x) --> sum[mj/rhoj*f(xj)*W(x-xj,h), j]
! grad f(x) --> -sum[mj/rhoj*f(xj)*grad W(x-xj,h), j]
!
! - time step
! (1) registration - set section (not implemented in this code)
! (2) set density
! (3) set pressure
! (4) Verlet step1 (t = t + dt/2)
! (5) set acceleration ax ay and power <-- x(t+dt/2),y(t+dt/2)
! (6) Verlet step2 (t = (t + dt/2) + dt/2)
! (7) set boundary
! goto (1)
!
MODULE sph2d
MODULE OPTION ARITHMETIC NATIVE
MODULE OPTION BASE 0
PUBLIC SUB setInitialCondition,timeEvolution,drawParticles
SHARE NUMERIC sysTime, dt, xMax, yMax, hh0, NNp
SHARE NUMERIC xx(5000),yy(5000) ! (xx(i),yy(i)) : position of i-th particle
SHARE NUMERIC vx(5000),vy(5000) ! (vx(i),vy(i)) : velocity of i-th particle
SHARE NUMERIC ax(5000),ay(5000) ! (ax(i),ay(i)) : total force/mass of i-th particle
SHARE NUMERIC hh(5000),mass(5000),cp(5000)
SHARE NUMERIC density(5000),pressure(5000),energy(5000),power(5000)
LET sysTime = 0.0 ! (s) system time
LET dt = 0.0001 ! (s) time step
LET xMax = 5.0 ! (m) x-Box size
LET yMax = 5.0 ! (m) y-Box size
LET hh0 = 0.25 ! (m) h of weight function
LET NNp = 400 ! number of smoothed particle
! ---------- set initial condition
EXTERNAL SUB setInitialCondition
! i,j, s
DECLARE EXTERNAL FUNCTION kernel,distance
RANDOMIZE
LET sysTime = 0.0
!set smoothed particles
FOR i=0 TO NNp-1
LET mass(i) = (0.029/0.0224)/4.0 !air
LET cp(i) = 1000.0 !air
LET xx(i) = 0.1*MOD(i,20)+0.5
LET yy(i) = 0.1*INT(i/20)+1.0
LET vx(i) = 0.0
LET vy(i) = 0.0
LET ax(i) = 0.0
LET ay(i) = 0.0
LET hh(i) = hh0
NEXT i
! set energy and power
FOR i=0 TO NNp-1
LET energy(i) = mass(i)*cp(i)*300 ! 300=temp(K)
LET power(i) = 0.0
NEXT i
!set density
FOR i=0 TO NNp-1
LET s = 0.0
FOR j=0 TO NNp-1
LET s = s + mass(j)*kernel(distance(i,j), (hh(i)+hh(j))/2.0)
NEXT j
LET density(i) = s
NEXT i
! set window
SET WINDOW 0,500, 0,500
END SUB
! ---------- time evolution
EXTERNAL SUB timeEvolution(nCalc)
! i
DECLARE EXTERNAL SUB timeStep
FOR i=0 TO NNp-1
CALL timeStep
LET sysTime = sysTime + dt
NEXT i
END SUB
EXTERNAL SUB timeStep
! i,j, dtv2,s,rr
DECLARE EXTERNAL SUB accCalc
DECLARE EXTERNAL FUNCTION kernel,distance
LET dtv2 = dt/2.0
!(2) set density
FOR i=0 TO NNp-1
LET s = 0.0
FOR j=0 TO NNp-1
LET s = s + mass(j)*kernel(distance(i,j), (hh(i)+hh(j))/2.0)
Next j
LET density(i) = s
NEXT i
!(3) SET pressure
FOR i=0 TO NNp-1
LET pressure(i) = 8.31*density(i)*(energy(i)/(mass(i)*cp(i)))
NEXT i
!(4) Verlet step1 (t = t + dt/2)
FOR i=0 TO NNp-1
LET vx(i) = vx(i) + dtv2*ax(i)
LET vy(i) = vy(i) + dtv2*ay(i)
LET xx(i) = xx(i) + vx(i)*dt
LET yy(i) = yy(i) + vy(i)*dt
NEXT i
!(5) set acceleration ax[] ay[] and power[]
CALL accCalc
!(6) Verlet step2 (t = t + dt/2 + dt/2)
FOR i=0 TO NNp-1
LET vx(i) = vx(i) + dtv2*ax(i)
LET vy(i) = vy(i) + dtv2*ay(i)
LET energy(i) = energy(i) + power(i)*dt
NEXT i
!(7) boundary
LET rr = 1.0
FOR i=0 TO NNp-1
if (xx(i) < 0.0) then
LET xx(i) = 0.0
LET vx(i) = -rr*vx(i)
LET vy(i) = rr*vy(i)
END IF
if (xx(i) > xMax) then
LET xx(i) = xMax
LET vx(i) = -rr*vx(i)
LET vy(i) = rr*vy(i)
end if
if (yy(i) < 0.0) then
LET yy(i) = 0.0
LET vx(i) = rr*vx(i)
LET vy(i) = -rr*vy(i)
END IF
if (yy(i) > yMax) then
LET yy(i) = yMax
LET vx(i) = rr*vx(i)
LET vy(i) = -rr*vy(i)
END IF
NEXT i
END SUB
! --- (5) set acceleration ax[] ay[] and power[]
EXTERNAL SUB accCalc
! i,j, ai,aj,b,rij,gradW,gradWx,gradWy
DECLARE EXTERNAL FUNCTION dwwvdr
FOR i=0 TO NNp-1
LET ax(i) = 0.0
LET ay(i) = 0.0
LET power(i) = 0.0
LET ai = pressure(i)/(density(i)*density(i))
FOR j=0 TO NNp-1
LET rij = SQR((xx(i)-xx(j))*(xx(i)-xx(j))+(yy(i)-yy(j))*(yy(i)-yy(j)))
IF (rij>0) THEN
LET aj = pressure(j)/(density(j)*density(j))
LET b = mass(j)*(ai+aj)
LET gradW = dwwvdr(rij, (hh(i)+hh(j))/2.0)
LET gradWx = gradW*(xx(i)-xx(j))/rij
LET gradWy = gradW*(yy(i)-yy(j))/rij
LET ax(i) = ax(i) - b*gradWx
LET ay(i) = ay(i) - b*gradWy
LET power(i) = power(i) + 0.5*b*((vx(i)-vx(j))*gradWx+(vy(i)-vy(j))*gradWy)
END IF
NEXT j
NEXT i
END SUB
! --- smoothed particle
EXTERNAL FUNCTION distance(i,j)
! x,y
LET x = xx(i)-xx(j)
LET y = yy(i)-yy(j)
LET distance = SQR(x*x+y*y)
END FUNCTION
EXTERNAL FUNCTION kernel(r,h)
! a,q,ret
LET ret = 0.0
LET q = r/h
LET a = 10.0/(7.0*3.141592)/(h*h)
if (q<1.0) then
LET ret = a*(1.0-1.5*q*q+0.75*q*q*q)
ELSEIF (q<2.0) THEN
LET ret = a*0.25*(2.0-q)*(2.0-q)*(2.0-q)
END IF
LET kernel = ret
END FUNCTION
EXTERNAL FUNCTION dwwvdr(r,h)
! a,q,ret
LET ret = 0.0
LET q = r/h
LET a = 10.0/(7.0*3.141592)/(h*h*h)
if (q<1.0) then
LET ret = a*(-3.0*q+2.25*q*q)
elseif (q<=2.0) then
LET ret = -a*0.75*(2.0-q)*(2.0-q)
END IF
LET dwwvdr = ret
END FUNCTION
! ---------- utility
! ---------- draw particles
EXTERNAL SUB drawParticles
DECLARE EXTERNAL FUNCTION systemTemprature
LET boxsize = 300
LET sc = boxsize/xMax
LET xp = 100
LET yp = 80
LET vScate = 100*dt
LET fScale = 1000*dt*dt/mass(0)
SET DRAW MODE HIDDEN
CLEAR
SET LINE COLOR 1 ! black : wall
PLOT LINES: xp,yp; xp+boxSize,yp; boxSize+xp,boxSize+yp; xp,boxSize+yp; xp,yp
!draw Balls, velocity and force
FOR i=1 TO nMolec
SET LINE COLOR 2 ! blue : molecule
DRAW circle WITH SCALE(sigma/2.0*sc)*SHIFT(xx(i)*sc+xp,yy(i)*sc+yp)
SET LINE COLOR 4 ! red : velocity
PLOT LINES: xx(i)*sc+xp,yy(i)*sc+yp;
PLOT LINES: (xx(i)+vx(i)*vScate)*sc+xp,(yy(i)+vy(i)*vScate)*sc+yp
!SET LINE COLOR 1 ! black : force
!PLOT LINES: xx(i)*sc+xp,yy(i)*sc+yp;
!PLOT LINES: (xx(i)+ffx(i)*fScale)*sc+xp,(yy(i)+ffy(i)*fScale)*sc+yp
NEXT i
SET TEXT HEIGHT 6
SET LINE COLOR 4 ! red : velocity
PLOT LINES: xp+0.6*boxsize,yp+boxsize+25;xp+0.6*boxsize+50,yp+boxsize+25
SET TEXT COLOR 4 ! red
PLOT TEXT, AT xp+0.6*boxsize+55,yp+boxsize+22: "velosity"
SET LINE COLOR 1 ! black : force
PLOT LINES: xp+0.6*boxsize,yp+boxsize+10;xp+0.6*boxsize+50,yp+boxsize+10
SET TEXT COLOR 1 ! black
PLOT TEXT, AT xp+0.6*boxsize+55,yp+boxsize+7: "force"
SET TEXT HEIGHT 10
SET TEXT COLOR 1 ! black
PLOT TEXT, AT 50, 50 ,USING "time =####.#### ":sysTime
PLOT TEXT, AT 50, 35 ,USING "N =####":NNp
PLOT TEXT, AT 50, 20 :"Ar in the box (2-dimensional molecular dynamics)"
SET TEXT HEIGHT 6
PLOT TEXT, AT xp,boxSize+2+yp ,USING "box size =##.# x ##.# (m)":xMax,yMax
SET DRAW MODE EXPLICIT
END SUB
END MODULE
|
|