cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cc This code uses the standard RK3, without modifying cc the u bar values ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc parameter (N=64, NN=2*N+1) implicit real*8(a-h,o-z) real x(0:NN), ubar(0:NN,0:3) real uapp1(0:NN), xL(0:NN,0:3), si(0:NN) complex abar(0:NN), al(0:NN) complex w, A, y, ss, ss1, ss2, F, sum, sum1 pi=4.0*atan(1.0) w=cmplx(0, 1) alpha=15.0*log(10.0) dx=2.0*pi/NN dx1=0.5*dx N2=0.0*N dt=0.005 ttime=0.0 Tmax=2.0 par=1.9 c compute the initial fourier coefficient do j=0, NN x(j)=j*dx if ( (x(j)-dx1) .lt. par .and. (x(j)+dx1) .gt. par) jp=j enddo do j=0, NN s1=cos(0.5*(x(j)+dx1)) s2=cos(0.5*(x(j)-dx1)) if (j .lt. jp) then ubar(j, 0)=-2.0*(s1-s2)/dx elseif (j .gt. jp) then ubar(j, 0)=2.0*(s1-s2)/dx else p1=-2.0*(cos(0.5*par)-s2)/dx p2=2.0*(s1-cos(0.5*par))/dx ubar(j, 0)=p1+p2 endif s1=-cos(x(j)+dx1) s2=cos(x(j)-dx1) c ubar(j, 0)=0.3+0.7*(s1+s2)/dx enddo 99 do NR=0, 2 do k=0, NN-1 k1=k-0.5*(NN-1) ss=0.0 do j=0, NN-1 p1=k1*x(j) s1=cos(p1) s2=-sin(p1) ss=ss+ubar(j,NR)*cmplx(s1,s2) enddo abar(k)=ss/NN enddo do k=0, NN-1 k1=k-0.5*(NN-1) if (k1 .eq. 0) then si(k)=1.0 else p1=0.5*k1*dx p2=sin(p1) si(k)=p2/p1 endif al(k)=abar(k)/si(k) enddo c compute the shock location N1=(1.0*N)**0.6 ss=(N1+1)*al(N+N1+1)/((N1+2)*al(N+N1+2)) y=-w*log(ss) ss=w*(N1+1)*y A=w*(N1+1)*al(N+N1+1)*exp(ss) if (NR .eq. 2) print *, 'haaaaaaaaaaa', A, y if (y1 .gt. 2.*pi .or. y1 .lt. 0.0) then y1=0.0 y=0.0 endif y1=real(y) y=real(y) A=real(A) do j=0, NN-1 xj=(j+0.5)*dx if (xj .le. y1 ) then F=-A*xj else F=A*(2*pi-xj) endif sum=0.0 do k=1, N if ( k .le. N2) then tau=1.0 else tau=exp(-alpha*(1.0*(k-N2)/N)**8) endif ps=k*y1 s1=cos(ps) s2=-sin(ps) ss1=cmplx(s1, s2) ps=k*xj s1=cos(ps) s2=sin(ps) ss2=cmplx(s1, s2) ss=(al(k+N)+w*A*ss1/k)*ss2 sum=sum+tau*ss k1=-k ps=k1*y1 s1=cos(ps) s2=-sin(ps) ss1=cmplx(s1, s2) ps=k1*xj s1=cos(ps) s2=sin(ps) ss2=cmplx(s1, s2) ss=(al(k1+N)+w*A*ss1/k1)*ss2 sum=sum+tau*ss enddo uapp1(j)=real(al(N)-A*(pi-y)+sum+F) enddo c Runge-Kutta do j=0, NN-1 ss1=uapp1(j) if (j .eq. 0) then ss2=uapp1(NN-1) else ss2=uapp1(j-1) endif xL(j, NR)=-0.5*(ss1*ss1-ss2*ss2)/dx if (NR .eq. 0) then ubar(j,1)=ubar(j,0)+dt*xL(j,0) endif if (NR .eq. 1) then sa=0.75*ubar(j,0)+0.25*ubar(j,1)+0.25*dt*xL(j,1) ubar(j,2)=sa endif if (NR .eq. 2) then sa=ubar(j,0)/3.0+ubar(j,2)*2/3.0+dt*xL(j,2)*2/3.0 ubar(j,3)=sa endif enddo enddo ttime=ttime+dt print *, ttime, ' time' if (ttime .lt. Tmax) then do j=0, NN ubar(j,0)=ubar(j, 3) enddo goto 99 endif do k=0, NN-1 k1=k-0.5*(NN-1) ss=0.0 do j=0, NN-1 p1=k1*x(j) s1=cos(p1) s2=-sin(p1) ss=ss+ubar(j,NR)*cmplx(s1,s2) enddo abar(k)=ss/NN al(k)=abar(k)/si(k) enddo do j=0, NN-1 xj=x(j) if (xj .le. y1 ) then F=-A*xj else F=A*(2*pi-xj) endif sum=0.0 do k=1, N if ( k .le. N2) then tau=1.0 else tau=exp(-alpha*(1.0*(k-N2)/N)**4) endif ps=k*y1 s1=cos(ps) s2=-sin(ps) ss1=cmplx(s1, s2) ps=k*xj s1=cos(ps) s2=sin(ps) ss2=cmplx(s1, s2) ss=(al(k+N)+w*A*ss1/k)*ss2 sum=sum+tau*ss k1=-k ps=k1*y1 s1=cos(ps) s2=-sin(ps) ss1=cmplx(s1, s2) ps=k1*xj s1=cos(ps) s2=sin(ps) ss2=cmplx(s1, s2) ss=(al(k1+N)+w*A*ss1/k1)*ss2 sum=sum+tau*ss enddo uapp1(j)=real(al(N)-A*(pi-y)+sum+F) enddo write (17, *) 'x=[' do j=0, NN-1 write (17, *) x(j) enddo write (17, *) '];' write (18, *) 'y=[' do j=0, NN-1 write (18, *) real(uapp1(j)) enddo write (18, *) '];' print *, A, y, ' A, y' stop end