====== Battle of the Winds: Inside the Code ======
Other than the indat file, we modified/added four files in VH1/src/Serial . These can be {{::battleofthewinds_code.zip|downloaded here}} and are explained on this page.
=== 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
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 [[group8|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
=== 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
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
=== 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 750x250 zones for our 12x4 AU problem.