subroutine evolve( umid, pmid ) ! Use umid and pmid from Riemann solver to update velocity, density, and total energy. ! Physical zones are from nmin to nmax. Zone boundary numbers run from nmin to nmax+1 !----------------------------------------------------------------------- ! GLOBALS use global use sweeps IMPLICIT NONE ! LOCALS INTEGER :: n REAL :: dtheta REAL, DIMENSION(maxsweep) :: umid, pmid, amid, uold, xa1, dvol1, upmid, dm, dtbdm REAL, DIMENSION(maxsweep) :: grav0, grav1, xa2, fict0, fict1, xa3 REAL, PARAMETER :: third = 1.0 / 3.0 !------------------------------------------------------------------------ ! grid position evolution do n = nmin-3, nmax + 4 dm (n) = r(n) * dvol(n) dtbdm(n) = dt / dm(n) xa1 (n) = xa(n) dvol1(n) = dvol(n) xa (n) = xa(n) + dt * umid(n) / radius upmid(n) = umid(n) * pmid(n) enddo xa1(nmin-4) = xa(nmin-4) xa1(nmax+5) = xa(nmax+5) do n = nmin-4, nmax+5 xa2(n) = xa1(n) + 0.5*dx(n) dx (n) = xa(n+1) - xa(n) xa3(n) = xa (n) + 0.5*dx(n) enddo ! Calculate forces using coordinates at t and at t+dt, note that ! fictitious forces at t+dt depend on updated velocity, but we ignore this call forces(xa2,grav0,fict0) call forces(xa3,grav1,fict1) ! Calculate dvolume and average area based on geometry of sweep select case (ngeom) case(0) do n = nmin-3, nmax+4 dvol(n) = dx(n) amid(n) = 1.0 enddo case(1) do n = nmin-3, nmax+4 dvol(n) = dx(n)*(xa(n)+.5*dx(n)) amid(n) = 0.5*(xa(n) + xa1(n)) enddo case (2) do n = nmin-3, nmax+4 dvol(n) = dx(n)*(xa(n)*(xa(n)+dx(n))+dx(n)*dx(n)*third) amid(n) = (xa(n)-xa1(n))*(third*(xa(n)-xa1(n))+xa1(n))+xa1(n)**2 enddo case (3) do n = nmin-3, nmax+4 dvol(n) = dx(n)*radius amid(n) = 1.0 enddo case (4) do n = nmin-3, nmax+4 dvol(n) = (cos(xa(n))-cos(xa(n+1)))*radius dtheta = xa(n) - xa1(n) if(dtheta .eq. 0.0) then amid(n) = sin(xa(n)) else amid(n) = (cos(xa1(n))-cos(xa(n)))/dtheta endif enddo case (5) do n = nmin-3, nmax+4 dvol(n) = dx(n) * radius amid(n) = 1.0 enddo case default write(*,*) 'Geometry', ngeom, ' not implemented.' stop end select do n = nmin-3, nmax+3 ! density evolution. lagrangian code, so all we have to do is watch the change in the geometry. r(n) = r(n) * ( dvol1(n) / dvol(n) ) r(n) = max(r(n),smallr) ! velocity evolution due to pressure acceleration and forces. uold (n) = u(n) u(n) = u(n) - dtbdm(n)*(pmid(n+1)-pmid(n))*0.5*(amid(n+1)+amid(n)) + 0.5*dt*(grav0(n)+fict0(n)+grav1(n)+fict1(n)) ! total energy evolution e(n) = e(n) - dtbdm(n)*(amid(n+1)*upmid(n+1) - amid(n)*upmid(n)) + 0.5*dt*(uold(n)*grav0(n) + u(n)*grav1(n)) q(n) = e(n) - 0.5*(u(n)**2+v(n)**2+w(n)**2) q(n) = max(q(n),smallp/(gamm*r(n))) enddo return end