Nuclear Talent
Course 7
NT4A

group8_code

# Battle of the Winds: Inside the Code

#### init.f90

We start by defining our computational grid, which is taken to be cylindrically symmetric:

ngeomx = 0 !cartesian = z
ngeomy = 1 !cylindrical = r
ngeomz = 0 !not used

nleftx = 1
nrightx= 1
nlefty = 0
nrighty= 1
nleftz = 0
nrightz= 0

xmin   = 0.0	! should be 0
xmax   = 12	! in AU
ymin   = 0.0	! should be 0
ymax   = 4	! in AU
zmin   = 0.0	! unused/2d problem
zmax   = 1.0	! unused/2d problem

Note that while VH-1 expects 3D cylindrical coordinates in the order r,phi,z, the order for a cylindrically symmetric 2D simulation is z,r.

Next, we define our input parameters:

!units:
![M]=10^12 g
![T]=yr
![L]=AU

gam    = 5. / 3.
gamm = gam - 1.0

! ambient material
pamb = 3.27e4 		!pressure
ramb = 5.60e3 		!density

! solar wind properties
zinj = 4		! z-coordinate of injection disk
widthinj = 1.6		! width of injection disk (radinj/widthinj <= r <= radinj)
vinj = 88.75 		! velocity of injected material
pinj = pamb 		! pressure of injected material
rinj = 20*4.54e4 	! density of injected material

! SNR shock
ncycleshock = 4500 	!number of cycles to wait for solar wind to settle before shock
vshock = 128.6 		!velocity of shock material
pshock = 1.23e8 	!pressure of shock material
rshock = 2.24e4 	!density of shock material

See our logbook for details on why we chose these units.

We then output all parameters to the hst file (code not shown here) and initialize our grid. As we need to stabilize the wind before introducing the shock, we only set up the solar wind source and the ambient medium here. As the solar wind is supposed to radiate spherically, we calculate the angle of the velocity of the solar wind for a given point on the injection radius.

do k = 1, kmax
do j = 1, jmax
do i = 1, imax
xphys=zxc(i) !float(i)/(xmax-xmin)
yphys=zyc(j) !float(j)/(ymax-ymin)

distcloud=sqrt((xphys-zinj)**2+(yphys)**2) !distance to cloud/wind center/'sun'
theta=atan2(yphys,xphys-zinj) !angle between z-axis and vector pointing from sun to current position

zro(i,j,k) = rinj
zpr(i,j,k) = pinj
zux(i,j,k) = vinj*cos(theta)
zuy(i,j,k) = vinj*sin(theta)
zro(i,j,k) = rinj
zpr(i,j,k) = pinj
zux(i,j,k) = 0
zuy(i,j,k) = 0
!ambient medium
else
zro(i,j,k) = ramb
zpr(i,j,k) = pamb
zux(i,j,k) = 0.
zuy(i,j,k) = 0.
endif
zuz(i,j,k) = 0.
zfl(i,j,k) = 0.
enddo
enddo
enddo

#### vhone.f90

As both the solar wind and the supernova are modeled as sources, we need to inject material in every timestep. This is done at the end of the main loop in vhone.f90.

  ! reset wind source each loop, i.e. keep it constant
do k = 1, kmax
do j = 1, jmax
do i = 1, imax
xphys=zxc(i) !float(i)/(xmax-xmin)
yphys=zyc(j) !float(j)/(ymax-ymin)

distcloud=sqrt((xphys-zinj)**2+(yphys)**2) !distance to cloud/wind center/'sun'
theta=atan2(yphys,xphys-zinj) !angle between z-axis and vector pointing from sun to current position

!stationary cloud
zpr(i,j,k) = pinj
zux(i,j,k) = vinj*cos(theta)
zuy(i,j,k) = vinj*sin(theta)
zro(i,j,k) = rinj
zro(i,j,k) = rinj
zpr(i,j,k) = pinj
zux(i,j,k) = 0
zuy(i,j,k) = 0
endif
enddo
enddo
enddo

!SNR coming in later
if(ncycle>ncycleshock) then
do j = 1, jmax
zpr(1,j,1) = pshock
zux(1,j,1) = vshock
zro(1,j,1) = rshock
enddo
endif

In addition, we had an issue with density accumulating along the symmetry axis. We think this is caused by issues with the symmetry/boundary conditions/ghost zones. Our fix (acknowledging Richard for the idea) is to set the densities in the innermost zones to the one in zone 3 in every timestep:

  !dirty fix: prevent the density from jumping at r=0 by setting the outermost gridpoint's density to the one of its radial neighbor
do i = 1, imax
zro(i,1,1) = zro(i,3,1)
zro(i,2,1) = zro(i,3,1)
enddo

#### solarwindmod.f90

As both init.f90 and vhone.f90 access the same variables containing information about the solar wind and the blast, we need these variables to be global. This is accomplished by introducing a new module called solarwind:

module solarwind
!=======================================================================
! global variables needed to initialize solar wind in init.f90 and keep it steady in vh90.f90
!=======================================================================
REAL :: zinj, radinj, vinj, pinj, rinj, widthinj, vshock, pshock, rshock
INTEGER :: ncycleshock

end module solarwind

#### zonemod.f90

As usual, the number of zones should be set according to the chosen geometry of the problem. If the number of zones is too small, the solar wind will be nowhere close to spherical. We ended up using 750×250 zones for our 12×4 AU problem.