#!/usr/bin/perl require "code.gplo"; #.... 1d heat equation, time dependent, Crank-Nicholson # make_f77(); system("./a.out "); # exit; #....... make plots, loop over output files $k=0; foreach $name (glob("tmpname-temp*")){ if($name=~/(\d{6})/){ $ist=$1; $pl1=" xlab=depth(m) lwidth=3 type=1 with=lines "; $ou='tmpx.pdf'; $ylo='set~logscale~y;set~format~y~"10^{%L}"'; $gpl=gpl_head("term=jpg font=Helvetica fsize=6 dx=0.5 title=1d~heat~equation~$ist"); $gpl=$gpl . gpl_xyplot("file=$name using=1:2 key=Temperature(C) yfix=0:100 ylab=Temperature $pl1 "); #$gpl=$gpl . gpl_xyplot("file=$name using=1:3 key=Conductivity ylab=kappa(W/m/K) lcol=rgb~'#dd9999' $pl1 "); #$gpl=$gpl . gpl_xyplot("file=$name using=1:4 key=Source ylab=Power(W/m^3) lcol=rgb~'#5555ff' yfix=1.0e-9:0.9 $pl1 more=$ylo "); gpl_final($gpl,"out=$ou"); ssy("convert -quality 100 -crop 800x400+0+540 tmp.jpg tmp_$k.jpg"); $k++; # ssy("xv tmp_$ist.jpg"); }} #....... make movie system("/bin/rm -fv movie.mp4"); #...... clean out stale files system("ffmpeg -i tmp_%d.jpg movie.mp4"); print STDERR "\n\n mpv --loop=inf movie.mp4\n\n"; sub ssy{ local($c)=@_; print STDERR "ssy:|$c|\n"; system($c); return; } sub make_f77{ $t=" c--------------------------------------------------- program heat1d c--------------------------------------------------- c time-dependent heat equation c JR 11/2017 c heat equation for vertical earth profile c x(i) is depth in meter c u(i) is temperature in C external xkappa,rho parameter(N=900) real x(N),T(N),Tnew(N),d(N),a(N),b(N),c(N) character*80 fname common /com1/ icase icase=2 if(icase.eq.2) then time0=0.0 yy=3600.0*24.0*365.0 timef=1000.0*yy x1=0.0 t1=10.0 !...... left boundary, earth surface xn=3000.0 tn=19.0 !...... right boundary, 3km depth, assuming 3k/km gradient dx=(xn-x1)/float(N-1) dt=7.5 !.... unstable dt=5.5 !.... stable dt=5.6 !.... unstable dt=100.0 !.... unstable dt=3600.0 dt=3600.0*24.0 !.... stable w/CR dt=10*3600.0*24.0 !.... stable w/CR cfl=0.5*dx*dx ck=xkappa(1000.0)/rhocv(1000.0) cfl=cfl/ck !..... adjust for real values write(0,*)'cfl,dt ',cfl/yy,dt dt=0.1*yy dtout=10.0*yy !...... initial values, step function do i=1,N x(i)=x1+float(i-1)*dx T(i)=t1+float(i-1)*(tn-t1)/float(N-1) !...... linear T Tnew(i)=T(i) enddo endif !....... test case if(icase.eq.1) then time0=0.0 timef=10.0 t1=-1.0 tn=1.0 x1=0.0 xn=10.0 dx=(xn-x1)/float(N-1) !...... initial values, step function do i=1,N x(i)=x1+float(i-1)*dx T(i)=t1 if(i.gt.N/2)T(i)=tn !..... step function Tnew(i)=T(i) enddo dt=0.1 dtout=0.5 endif c...... time step loop time=time0 kstep=0 outnext=time0 100 continue c..... set up equations do i=2,N-1 xx=x(i) xm=xx-0.5*dx xp=xx+0.5*dx rcv=rhocv(xx) r=dt/dx/dx/rcv/2.0 xkm=r*xkappa(xm) xkp=r*xkappa(xp) xk=xkm+xkp !...... right hand side and the 3 matrix diagonals d(i)=(1.0-xk)*T(i)+xkp*T(i+1)+xkm*T(i-1) +dt*qsource(xx,time) !... RHS a(i)=-xkm !.... lower diagonal b(i)=1.0+xk ! diagonal c(i)=-xkp ! upper diagonal enddo !...... adjust for boundary conditions d(2)=d(2) + t1*xkappa(x1+0.5*dx)*r d(N-1)=d(N-1)+ tn*xkappa(xn-0.5*dx)*r ! do i=2,N-1 ! write(0,*)i,a(i),b(i),c(i),d(i),r ! enddo !...... solve tridiagonal system ! Tnew(1) and Tnew(N) are already given by the boundary conditions ! Note we do the fortran 'pass by pointer' trick here call thomas(N-2,a(2),b(2),c(2),d(2),Tnew(2)) Tnew(1)=T(1) Tnew(N)=T(N) !..... not really needed, just for clarity ! do i=1,N ! write(0,*)i,T(i),Tnew(i) ! enddo !...... copy back do i=1,N T(i)=Tnew(i) enddo time=time+dt !........ output some points to look at progress write(0,*)time/yy,T(300),T(440),T(455),T(600) !...... output every so often if(time.gt.outnext) then k=int(dtout+0.5); m=int(time); m=m/k; iout=m*k; write(0,*)' m,k,iout ',m,k,iout fname=' ' write(fname,'(a,i6.6)')'tmpname-temp.',int(time/yy) call wout(n,x,T,fname) outnext=time+dtout endif !...... next step if still time left kstep=kstep+1 if(kstep.gt.10000)stop if(time.lt.timef) goto 100 stop end c--------------------------------------------------- subroutine wout(n,x,y,fname) c--------------------------------------------------- real x(n),y(n) character*(*) fname write(0,*)'writing ... ',TRIM(fname) open(10,file=fname) do i=1,n write(10,*)x(i),y(i) enddo close(10) return end c--------------------------------------------------- subroutine thomas(n,a,b,c,d,x) c--------------------------------------------------- real a(n),b(n),c(n),d(n),x(n) !... note we overwrite c and d common /com1/ icase c(1)=c(1)/b(1) d(1)=d(1)/b(1) do i=2,n z=1.0/(b(i)-a(i)*c(i-1)) c(i)=c(i)*z d(i)=(d(i)-a(i)*d(i-1))*z enddo x(n)=d(n) do i=n-1,1,-1 x(i)=d(i)-c(i)*x(i+1) enddo return end c--------------------------------------------------- real function rhocv(x) c--------------------------------------------------- c returns material rho*Cv common /com1/ icase cv=2000 !..... heat capacity in J/K/kg rho=2750 !..... density of granite in kg/m^3 rhocv=rho*cv !..... in J/K/m^3 if(icase.eq.1)rhocv=1.0 !...... for testing and stbility demo return end c--------------------------------------------------- real function xkappa(x) c--------------------------------------------------- c returns thermal conductivity in W/m/K common /com1/ icase xkappa=1.8 !... granite if(x.lt.1000.0) xkappa=2.2 !.... sandstone if(x.lt.500.0) xkappa=1.0 !.... limestone if(x.lt.150.0) xkappa=0.5 !.... soil/gravel if(icase.eq.1)xkappa=1.0 !.... testing and demo return end c--------------------------------------------------- real function qsource(x,time) c--------------------------------------------------- c..... returns heating rate in W/m^3 common /com1/ icase yy=3600.0*24.0*365.0 qsource=0.0 ! qsource=1.0e-8 ! if(x.gt.1000.0)qsource=0.5e-6 !...... granite if((x.gt.800.0).and.(x.lt.825.0)) then qsource=3.0e-7*exp(-time/yy/30.0) !...... diluted radioactive waste after 10s of years ! write(0,*) 'qsource: ',x,qsource endif if(icase.eq.1)qsource=0.0 return end "; system("/bin/rm -f tmp.f a.out"); #...... clean up open(FF,">tmp.f"); print FF $t; close(FF); #...... write code to temporary file system("gfortran tmp.f"); if( -e "a.out" ){ system(" size a.out "); return; } else{ print STDERR "compilation failed!\n"; exit; } }