Boundary conditions

The computational domain

We used 4 ghost cells at each boundary of the domain. The computational domain is thus defined by

• nx+8 cells in the r-direction, and ny+8 cells in the theta-direction

The ghost cells

At each time step, we compute the solution of equations inside the hole computational domain. After what we adjust variables in the ghost cells.

• One assume that r(0) is the radius of the first ghost cell. Thus the index of the ghost cells are i = 0, 1, 2, 3, nx+5, nx+6, nx+7, nx+8.
• In these ghost cells in the r-direction, we compute the average value of ug, vg, pg, rosg, over the ny azimuthal cells. One call them ug_ave(i), vg_ave(i), pg_ave(i), rosg_ave(i)
• We actualize variables in all the ghost cells

do i = 0, 1, 2, 3
   do j = 0, ny+8

 !--- first cells ----
     ug(i,j) = 2**(-4+i)*(ug(i,j) - ug_ave(i))
     vg(i,j) = 2**(-4+i)*(vg(i,j) - vg_ave(i)) + dsqrt(1.D0/r(i) + r(i)**(3/2) * (-2)*pref*r(i)**(-2)) ! this is the keplerian velocity vk^2 = 1/r + r/rosg d/dr p
     pg(i,j) = 2**(-4+i)*(pg(i,j) - pg_ave(i)) + pref*r(i)**(-2)  ! pref = 1.125D-3
     rosg(i,j)= 2**(-4+i)*(rosg(i,j)- rosg_ave(i)) + r(i)**(-3/2) 

 !--- last cells ---
     ug(nx+8-i,j) = 2**(-4+i)*(ug(nx+8-i,j) - ug_ave(nx+8-i))
     vg(nx+8-i,j) = 2**(-4+i)*(vg(nx+8-i,j) - vg_ave(nx+8-i)) + dsqrt(1.D0/r(nx+8-i) + r(nx+8-i)**(3/2) * (-2)*pref*r(nx+8-i)**(-2))
     pg(nx+8-i,j) = 2**(-4+i)*(pg(nx+8-i,j) - pg_ave(nx+8-i)) + pref*r(nx+8)**(-2)
     rosg(nx+8-i,j)= 2**(-4+i)*(rosg(nx+8-i,j)- rosg_ave(nx+8-i)) + r(nx+8-i)**(-3/2) 
    
   end
end

This permits to absorb density waves, and to keep the power low in density and pressure.

• Finally we use periodical conditions in the theta direction.