Other than the indat file, we modified/added four files in VH1/src/Serial . These can be downloaded here and are explained on this page.

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
radinj = 0.25		! 1e18 ! injection disk radius
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
 
    !injection radius
    if(radinj/widthinj < distcloud .and. distcloud < radinj .and. j>0) then
      zro(i,j,k) = rinj
      zpr(i,j,k) = pinj
      zux(i,j,k) = vinj*cos(theta)
      zuy(i,j,k) = vinj*sin(theta)
    !fixed density inside injection radius
    elseif(distcloud <= radinj/widthinj) then
      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

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
	if(radinj/widthinj < distcloud .and. distcloud < radinj .and. j>0) then
	  zpr(i,j,k) = pinj
	  zux(i,j,k) = vinj*cos(theta)
	  zuy(i,j,k) = vinj*sin(theta)
	  zro(i,j,k) = rinj
	elseif(distcloud <= radinj/widthinj) then
	  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

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

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.