group8_code

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.

group8_code.txt · Last modified: 2014/06/06 11:15 by bartl