solution.py

###############################################################################
 # This file is part of SWIFT.
 # Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
 #                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
 # 
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU Lesser General Public License as published
 # by the Free Software Foundation, either version 3 of the License, or
 # (at your option) any later version.
 # 
 # This program is distributed in the hope that it will be useful,
 # but WITHOUT ANY WARRANTY; without even the implied warranty of
 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 # GNU General Public License for more details.
 # 
 # You should have received a copy of the GNU Lesser General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 # 
 ##############################################################################

import random
from numpy import *

# Computes the analytical solution of the Gresho-Chan vortex
# The script works for a given initial box and background pressure and computes the solution for any time t (The solution is constant over time).
# The code writes five files rho.dat, P.dat, v.dat, u.dat and s.dat with the density, pressure, internal energy and
# entropic function on N radial points between r=0 and r=R_max.

# Parameters
rho0 = 1.         # Background Density
P0 = 0.           # Background Pressure
gamma = 5./3.     # Gas polytropic index
N = 1000          # Number of radial points
R_max = 1.        # Maximal radius

# ---------------------------------------------------------------
# Don't touch anything after this.
# ---------------------------------------------------------------

r = arange(0, R_max, R_max / N)
rho = ones(N)
P = zeros(N)
v = zeros(N)
u = zeros(N)
s = zeros(N)


for i in range(N):
    if r[i] < 0.2:
        P[i] = P0 + 5. + 12.5*r[i]**2
        v[i] = 5.*r[i]
    elif r[i] < 0.4:
        P[i] = P0 + 9. + 12.5*r[i]**2 - 20.*r[i] + 4.*log(r[i]/0.2)
        v[i] = 2. -5.*r[i]
    else:
        P[i] = P0 + 3. + 4.*log(2.)
        v[i] = 0.
    rho[i] = rho0
    s[i] = P[i] / rho[i]**gamma
    u[i] = P[i] /((gamma - 1.)*rho[i])

savetxt("rho.dat", column_stack((r, rho)))
savetxt("P.dat", column_stack((r, P)))
savetxt("v.dat", column_stack((r, v)))
savetxt("u.dat", column_stack((r, u)))
savetxt("s.dat", column_stack((r, s)))