( ) 1
1 # ls scalar Makefile anime.pro main.f pldt.pro pldtps.pro rddt.pro Fortran Fortran main.f Fortran 1.1 (make) scalar make main.o a.out out.dat # cd scalar # make f77 -c -o main.o main.f main.f: MAIN: f77 -o a.out main.o./a.out write step= 0 time= 0.000E+00 write step= 50 time= 0.125E+02 write step= 100 time= 0.250E+02 ### normal stop ### # ls Makefile anime.pro main.o pldt.pro rddt.pro a.out* main.f out.dat pldtps.pro 1.2 (out.dat) out.dat 1 (jx) (nx) 2 2
time step (ns) (time) 3 102 x (x) (u) 103 Fortran main.f 53, 55, 59 Format # head out.dat 100, 3 0, 0.00 1.0, 1.0000000 2.0, 1.0000000 3.0, 1.0000000 ( ) 1.3 IDL IDL ( idl 1.3.1 IDL (idl) idl # idl IDL IDL Version 5.5 (sunos sparc). (c) 2001, Research Systems, Inc. Installation number: XXXXX. Licensed for use by: XXXXX IDL> 1.3.2 (.r rddt) rddt.pro out.dat idl IDL>.r rddt.r run 3
1.3.3 (.r pldt) pldt.pro IDL>.r pldt 1: scalar 1.3.4 IDL (exit) exit IDL IDL> exit 4
2 2.1 FTCS 72 88 c---------------------------------------------------------------------- c solve equation c c ftcs - start >>> do j=1,jx-1 f(j)=0.5*cs*(u(j+1)+u(j)) enddo f(jx)=f(jx-1) do j=2,jx-1 u(j)=u(j)-dt/dx*(f(j)-f(j-1)) enddo u(1)=u(2) u(jx)=u(jx-1) c ftcs - end >>> 2.2 (jx) 5 parameter jx parameter (jx=100) 5
2.3 (nstop, nskip) 14 nstop 15 nskip c time control parameters nstop=100 nskip = 50 2.4 CFL (safety) CFL 68 (safety) nskip safety c obtain time spacing safety=0.25 3 (.r anime) anime.pro IDL IDL>.r anime nskip 1 anime.pro window anime.pro 6
main.f c====================================================================== c array definitions c====================================================================== implicit real*8 (a-h,o-z) parameter (jx=100) dimension x(1:jx),u(1:jx),f(1:jx) c====================================================================== c prologue c====================================================================== c time control parameters nstop=100 nskip = 50 c---------------------------------------------------------------------- c initialize counters time = 0.0 ns = 0 nx = nstop/nskip+1 c---------------------------------------------------------------------- c Set initial condition c---------------------------------------------------------------------- pi=4.*atan(1.0) c grid dx=1.0 x(1)=dx do j=1,jx-1 x(j+1)=x(j)+dx enddo 7
c c c c variable do j=1,jx/2 u(j)= 1.0 enddo do j=jx/2+1,jx u(j)= 0.0 enddo velocity cs=1.0 c----------------------------------------------------------------------- c Output initial condition c write(6,103) ns,time 103 format (1x, write, step=,i8, time=,e10.3) open(unit=10,file= out.dat,form= formatted ) write(10,100) jx,nx 100 format(i5,,,i5) write(10,101) ns,time 101 format (i5,,,f6.2) do j=1,jx write(10,102) x(j),u(j) enddo 102 format(f5.1,,,f10.7) c====================================================================== c time integration c====================================================================== 1000 continue ns = ns+1 c---------------------------------------------------------------------- c obtain time spacing safety=0.25 dt=safety*dx/cs time=time+dt 8
c---------------------------------------------------------------------- c solve equation c c ftcs - start >>> do j=1,jx-1 f(j)=0.5*cs*(u(j+1)+u(j)) enddo f(jx)=f(jx-1) do j=2,jx-1 u(j)=u(j)-dt/dx*(f(j)-f(j-1)) enddo u(1)=u(2) u(jx)=u(jx-1) c ftcs - end >>> c---------------------------------------------------------------------- c data output if (mod(ns,nskip).eq.0) then write(6,103) ns,time write(10,101) ns,time do j=1,jx write(10,102) x(j),u(j) enddo endif if (ns.lt. nstop) goto 1000 close(10) *======================================================================= write(6,*) ### normal stop ### end 9
rddt.pro ; rddt.pro openr,1, out.dat readf,1,jx,nx ; define array ns=intarr(nx) t=fltarr(nx) x=fltarr(jx) u=fltarr(jx,nx) ; temporary variables for read data ns_and_t=fltarr(2,1) x_and_u=fltarr(2,jx) for n=0,nx-1 do begin readf,1,ns_and_t readf,1,x_and_u ns(n)=fix(ns_and_t(0,0)) t(n)=ns_and_t(1,0) u(*,n)=x_and_u(1,*) endfor close,1 free_lun,1 x(*)=x_and_u(0,*) delvar,ns_and_t,x_and_u help end 10
pldt.pro!x.style=1!y.style=1!p.charsize=1.4 plot,x,u(*,0),xtitle= x,ytitle= u,linest=1,yrange=[-1,3],xrange=[0,100] for n=1,nx-1 do begin oplot,x,u(*,n) oplot,x,u(*,n),psym=4 endfor end anime.pro!x.style=1!y.style=1!p.charsize=1.4 window,xsize=480,ysize=480 xinteranimate,set=[480,480,nx] for n=0,nx-1 do begin plot,x,u(*,n),xtitle= x,ytitle= u,yrange=[-1,3],xrange=[0,100] oplot,x,u(*,n),psym=4 xinteranimate,frame=n,window=0 endfor xinteranimate end 11
1 1 (p273 p282) FTCS j = 1,...50 u j = 1 j = 51,...100 u j = 0 ν = c t/ x = 0.25 1. Lax-Wendroff 2. 1 3. minmod ( 1.12, 1.10) 1 u t + c u x = 0 (1) u n+1 j = u n j t x (f j+1/2 n fj 1/2) n (2) FTCS f n j+1/2 = 1 2 (f j+1 + f j ) = 1 2 c(u j+1 + u j ) (3) p21 p276 3 3 2 2 u 1 u 1 0 0-1 0 20 40 60 80 100 x -1 0 20 40 60 80 100 x 2: Lax-Wendroff + 1 12
2 Burgers 1 Burgers u t + x ( ) u 2 = 0 (4) 2 1 ( 1.15, 1.16) p23 f n j+1/2 = 1 2 {( uj+1 2 2 + u j 2 ) 1 } 2 2 u j+1 + u j (u j+1 u j ) (5) 3 1 1 1 u t = u κ 2 (6) x 2 FTCS cs kappa c c c variable do j=1,jx u(j)= exp(-(((x(j)-x(jx/2))/5.)**2)) enddo kappa kappa=1.0 13
14
1 CANS Coordinated Astronomical Numerical Software) 2 CANS CANS CANS CANS 1 3 3 CANS CANS CANS 4 CANS CANS Lax-Wendroff 15
Roe CIP-MOCCT CANS 1 2 3 MPI CANS/MPI CANS 1 CANS1D 2 3 CANS2D, CANS3D 5 CANS Makefile Makefile cans1d/ CANS 1D (1 ) Fortran cans2d/ CANS 2D (2 ) Fortran cans3d/ CANS 3D (3 ) Fortran cansnc/ idl/ IDL cans-1.0 cans3d 6 CANS Web http://www.astro.phys.s.chiba-u.ac.jp/netlab/ CANS http://www.astro.phys.s.chiba-u.ac.jp/netlab/astro/index.html 16
1 CANS 1D CANS 1D Lax-Wendroff hdmlw # ls hdmlw Makefile mlw_ht.f mlw_m3_g.f mlw_m_g.f mlwfull.f README mlw_ht_c.f mlw_m3t.f mlw_mt.f mlwhalf.f Readme.tex mlw_ht_cg.f mlw_m3t_c.f mlw_mt_c.f mlwsrcf.f mlw_a.f mlw_ht_g.f mlw_m3t_cg.f mlw_mt_cg.f mlwsrch.f mlw_h.f mlw_m.f mlw_m3t_g.f mlw_mt_cgr.f mlw_h_c.f mlw_m3.f mlw_m_c.f mlw_mt_g.f mlw_h_cg.f mlw_m3_c.f mlw_m_cg.f mlw_rh.f mlw_h_g.f mlw_m3_cg.f mlw_m_cgr.f mlwartv.f ( ) md # ls md_shktb Makefile bnd.f pldt.pro README cipbnd.f rddt.pro Readme.pdf main.f shktb_analytic.pro Readme.tex main.pro anime.pro model.f Fortran Fortran.f Fortran README Readme.pdf HTML htdocs web browser Web 17
1.1 (make) cans1d cansnc 2 make cans1d make libcans1d.a cansnc make libcansnc.a # cd cans1d # make cd hdmlw; make f77 -O -c mlwfull.f ( ) touch.update # cd../cansnc # make ( ) cans1d cansnc make cans1d cans2d cans3d cansnc 18
1.2 (make) make (md shktb) params.txt **.dac # cd md_shktb # make f77 -O -c main.f f77 -O -c model.f f77 -O -c bnd.f f77 -O -c cipbnd.f f77 -o a.out main.o model.o bnd.o cipbnd.o -L.. -lcans1d./a.out write step= 0 time= 0.000E+00 nd = 1 write step= 75 time= 0.505E-01 nd = 2 write step= 154 time= 0.100E+00 nd = 3 write step= 221 time= 0.142E+00 nd = 4 stop step= 221 time= 0.142E+00 ### normal stop ### 1.3 IDL IDL ( idl 1.3.1 IDL (idl) idl # idl IDL IDL Version 5.5 (sunos sparc). (c) 2001, Research Systems, Inc. Installation number: XXXXX. Licensed for use by: XXXXX IDL> 19
1.3.2 (.r rddt) rddt.pro.r run IDL>.r rddt help IDL> help ( ) GM FLOAT = 1.40000 IX LONG = 208 NX LONG = 4 PR FLOAT = Array[208, 4] PR1 FLOAT = 0.100000 RO FLOAT = Array[208, 4] RO1 FLOAT = 0.125000 T FLOAT = Array[4] TE FLOAT = Array[208, 4] VX FLOAT = Array[208, 4] X FLOAT = Array[208] Compiled Procedures: $MAIN$ NCGETO1 NCGETOS NCGETS1 NCGETSS ( ) PR( ) RO( ) TE( ) VX( ) 2 ( ) X( ) T( ) 1 help idl PR pr 1.3.3 (.r pldt) pldt.pro IDL>.r pldt 1.3.4 (.r anime) idl anime.pro anime.pro 20
3: md shktb IDL>.r anime 1.3.5 IDL (exit) exit IDL 7 IDL> exit 21
1 1. (md itshktb) IDL rddt.pro pldt.pro 2. (md shktb) 3. (md shkform) anime.pro 4. MHD (md mhdshktb) 4: md mhdshktb Fortran make params.txt ***.dac make clean 22
2 CANS 1D 3 main.f, model.f, bnd.f md *** # cp -r md_shktb md_shktb1 2.1 (main.f) main.f 2.1.1 (ix) (md shktb) (ix) main.f # cd md_shktb # ls Makefile bnd.f main.o pldt.pro Makefile-nc Readme.pdf bnd.o model.f shktb_analytic.pro Makefile-pg Readme.ps cipbnd.f model.o pr.dac Makefile-pgnc cipbnd.o params.txt rddt.pro t.dac README a.out* rdnc.pro vx.dac anime.pro main.f pldt.f ro.dac x.dac main.f (1 5 ) c====================================================================== c array definitions c====================================================================== implicit real*8 (a-h,o-z) parameter (ix=1026) 1026 2 407 make 2 1 4 2 8 23
2.1.2 (tend, dtout) (tend) (dtout) main.f (26 32 ) c---------------------------------------------------------------------- c time control parameters c nstop : number of total time steps for the run dtout=0.05 tend=0.14154 nstop=1000000 0.01 (md shktb) 24
2.2 (model.f) model.f model.f (ro) (vx, vy) (pr) (bx, by) 2.2.1 γ (gm) Sedov γ (gm) model.f model.f (13 16 ) c---------------------------------------------------------------------- c parameters c---------------------------------------------------------------------- gm=5./3. 5./3. 7./5. gm=7./5. make idl :Sedov (md sedov) 25
2.3 (main.f) 2.3.1 Roe MHD (md mhdshktb) Roe main.f (120 ) c solve hydrodynamic equations c hdmlw - start >>> qav=3.d0 call mlw_m(ro,pr,vx,vy,by,bx,bxm,dt,qav,gm,dx,dxm,ix) c hdmlw - end <<< c hdroe - start >>> c call roe_m(ro,pr,vx,vy,by,bx,bxm,dt,gm,dx,ix) c hdroe - end <<< c hdcip - start <<< c cvis=3.00d0 c call cip_m(ro,pr,vx,vy,by,te,vxm,rodx,tedx,vxdxm,vydx c &,bx,bxm,dt,cvis,gm,dx,dxm,ix) c hdcip - end <<< mlw m Lax-Wendroff MHD roe m Lax-Wendroff MHD call mlw m call roe m 26
c solve hydrodynamic equations c hdmlw - start >>> c qav=3.d0 c call mlw_m(ro,pr,vx,vy,by,bx,bxm,dt,qav,gm,dx,dxm,ix) c hdmlw - end <<< c hdroe - start >>> call roe_m(ro,pr,vx,vy,by,bx,bxm,dt,gm,dx,ix) c hdroe - end <<< c hdcip - start <<< c cvis=3.00d0 c call cip_m(ro,pr,vx,vy,by,te,vxm,rodx,tedx,vxdxm,vydx c &,bx,bxm,dt,cvis,gm,dx,dxm,ix) c hdcip - end <<< Roe Lax-Wendroff safety main.f 2.3.2 CIP MHD (md mhdshktb) CIP Roe main.f CIP 3 CIP hdcip c hdcip - start >>> dimension te(ix),vxm(ix),rodx(ix),tedx(ix),vxdxm(ix),vydx(ix) c hdcip - end <<< c hdcip - start >>> call ciprdy_m(te,vxm,rodx,tedx,vxdxm,vydx,ro,pr,vx,vy &,gm,dx,dxm,ix) call cipbnd(margin,ro,te,vxm,vy,by,rodx,tedx,vxdxm,vydx,ix) c hdcip - end <<< 27
c solve hydrodynamic equations c hdmlw - start >>> c qav=3.d0 c call mlw_m(ro,pr,vx,vy,by,bx,bxm,dt,qav,gm,dx,dxm,ix) c hdmlw - end <<< c hdroe - start >>> c call roe_m(ro,pr,vx,vy,by,bx,bxm,dt,gm,dx,ix) c hdroe - end <<< c hdcip - start <<< cvis=3.00d0 call cip_m(ro,pr,vx,vy,by,te,vxm,rodx,tedx,vxdxm,vydx &,bx,bxm,dt,cvis,gm,dx,dxm,ix) c hdcip - end <<< Roe Lax-Wendroff safety main.f 28
2.4 (bnd.f) bnd.f md shktb bnd.f call bdsppx(0,margin,ro,ix) call bdsppx(0,margin,pr,ix) call bdspnx(0,margin,vx,ix) call bdsppx(1,margin,ro,ix) call bdsppx(1,margin,pr,ix) call bdspnx(1,margin,vx,ix) bdsppx bdspnx bdsppx bdspnx bdsppx(,,, ) 4 3 3 margin margin Lax-Wendroff 1 Roe 2 CIP 4 3 (ro) (pr) (vx) bdspnx bdperx 4 bc/readme 29
CANS 1D (cans-1.0 ) hdmlw/ hdroe/ hdcip/ MHD Lax-Wendroff Roe CIP bc/ common/ cndsor/ 1 Red Black SOR cndbicg/ 1 BICG htcl/ selfg/ md_***/ README md_***/readme md_***/readme.pdf PDF md_***/readme.ps PS 30
md_advect/ [ ] md_shktb/ [ ] md_shkform/ [ ] md_strongshk/ [ ] md_itshktb/ [ ] md_mhdshktb/ MHD [MHD] md_itmhdshktb/ MHD [ MHD] md_awdecay/ Alfven [ 3 MHD] md_sedov/ Sedov [ ] md_thinst/ [ ] md_cndtb/ [ ] md_cndsp/ [ ] md_spicule/ [MHD ] md_wind/ Parker [ ] md_mhdwind/ MHD Weber-Davis [MHD ] md_diskjet/ MHD [MHD ] md_flare/ [ ] md_cloud/ [ ] md_clsp/ [ ] md_relshk/ [ ] bdcnsx qq(1)=q0 bdfrex qq(1)=qq(2) bdfrdx qq(1)=qq(2)-dqq(2)*dx bdperx qq(1)=qq(ix) bdsppx qq(1)= qq(2) bdspnx qq(1)=-qq(2) bdsmpx qq(1)= qq(3) bdsmnx qq(1)=-qq(3) 31
1 CANS 2D, 3D CANS 2D CANS 3D CANS 1D MPI Fortran mdp CANS 1D 1.1 (make) cans2d cansnc make cans2d make libcans2d.a # cd cans2d # make cd hdmlw; make f77 -c mlwfull.f mlwfull.f: mlwfull: f77 -c mlwhalf.f mlwhalf.f: mlwhalf: ( ) CANS 1D cansnc make libcansnc.a ( CANS 1D ) # cd../cansnc # make ( ) 32
CANS 3D cans3d make libcans3d.a # cd cans3d # make cd hdmlw; make f77 -c mlwfull.f mlwfull.f: mlwfull: f77 -c mlwhalf.f mlwhalf.f: mlwhalf: ( ) CANS 1D cans1d cans2d cans3d cansnc make 33
1.2 (make) make Kelvin-Helmholtz (md kh) # cd md_kh # make f77 -c main.f main.f: MAIN: f77 -c model.f model.f: model: f77 -c bnd.f bnd.f: bnd: f77 -o a.out main.o model.o bnd.o \ -L../.. -lcans2d -lcansnc./a.out write step= 0 time= 0.000E+00 nd = 1 write step= 83 time= 0.100E+01 nd = 2 write step= 167 time= 0.201E+01 nd = 3 write step= 251 time= 0.301E+01 nd = 4 write step= 336 time= 0.401E+01 nd = 5 write step= 421 time= 0.500E+01 nd = 6 write step= 510 time= 0.600E+01 nd = 7 write step= 605 time= 0.701E+01 nd = 8 write step= 707 time= 0.801E+01 nd = 9 write step= 817 time= 0.900E+01 nd = 10 write step= 940 time= 0.100E+02 nd = 11 stop step= 940 time= 0.100E+02 ### normal stop ### CANS 2D CANS 1D 34
1.3 IDL 1.3.1 IDL (idl) idl # idl 1.3.2 (.r rddt) IDL>.r rddt 1.3.3 (.r pldt) IDL>.r pldt % Compiled module: $MAIN$. Plot columns \& rows? : 2,2 Variable for color-maps? (ro,pr,te) : ro start step? : 0 2 2 2 ro pr te (vx, vy, vz) (bx, by, bz) 0 1.3.4 (.r anime) IDL>.r anime 35
1.3.5 (loadct,39 xloadct) loadct,39 pldt.pro IDL> loadct,39 loadct 0 40 0 B-W Linear( ) 39 rainbow+white( + ) xloadct IDL> xloadct 1.3.6 IDL (exit) IDL IDL> exit 36
CANS 2D (cans-1.0 ) hdmlw/ MHD Lax-Wendroff bc/ common/ cndsor/ 1 Red Black SOR cndbicg/ 1 BICG htcl/ selfgmg/ Multigrid Iteration commonmpi/ (MPI) cndsormpi/ 1 Red Black SOR (MPI) md_***/ mdp_***/ [ ] [ ] ( ) md shktb mdp shktb 37
md_advect/ [ ] md_shktb/ [ ] md_itshktb/ [ ] md_mhdshktb/ MHD [MHD] md_itmhdshktb/ MHD [ MHD] md_mhd3shktb/ 3 MHD [3 MHD] md_awdecay/ Alfven [ 3 MHD] md_thinst/ [ ] md_cndtb/ [ ] md_mhdcndtb/ MHD [MHD ] md_shkref/ [ ] md_kh/ Kelvin-Helmholtz [ ] md_rt/ Rayleigh-Taylor [ ] md_sndwave/ [ ] md_mhdwave/ MHD [MHD] md_mhdkh/ MHD Kelvin-Helmholtz [MHD] md_mhd3kh/ 3 MHD Kelvin-Helmholtz [3 MHD] md_recon/ [MHD ] md_recon3/ 3 [3 MHD ] md_efr/ Parker [MHD ] md_corjet/ [MHD ] md_mri/ [MHD Coriolis ] md_mricyl/ [MHD Coriolis ] md_reccnd/ [MHD ] md_cndsp/ [ ] md_sedov/ Sedov [ ] md_jetprop/ [ ] md_mhdsn/ MHD [MHD ] md_diskjet/ [MHD ] md_cme/ Low [MHD ] md_cloud/ [ ] md_mhdcloud/ MHD [ MHD ] md_mhdgwave/ MHD [MHD] md_parker/ Parker [MHD ] Version 3.2 (2003/08/08) 38
. Craig-Henton. ( )... MHD Craig-Henton Craig& Henton (1995) 2002 Craig-Henton Hirose et al. 2004, ApJ, 610,1107 ( ) MHD MHD CIP 2 MHD
MHD MPI MHD SUN Fujitsu VPP5000 MHD IMF VRML 3D-MHD, Modified Leap-Frog scheme, MPI