|
| 1 | +from __future__ import print_function |
| 2 | + |
| 3 | +import sys |
| 4 | +import numpy |
| 5 | +import mesh.patch as patch |
| 6 | +from util import msg |
| 7 | + |
| 8 | + |
| 9 | +def init_data(my_data, rp): |
| 10 | + """ initialize the Gresho vortex problem """ |
| 11 | + |
| 12 | + msg.bold("initializing the Gresho vortex problem...") |
| 13 | + |
| 14 | + # make sure that we are passed a valid patch object |
| 15 | + if not isinstance(my_data, patch.CellCenterData2d): |
| 16 | + print("ErrrrOrr: patch invalid in bubble.py") |
| 17 | + print(my_data.__class__) |
| 18 | + sys.exit() |
| 19 | + |
| 20 | + # get the density and velocities |
| 21 | + dens = my_data.get_var("density") |
| 22 | + xmom = my_data.get_var("x-momentum") |
| 23 | + ymom = my_data.get_var("y-momentum") |
| 24 | + ener = my_data.get_var("energy") |
| 25 | + |
| 26 | + # grav = rp.get_param("compressible.grav") |
| 27 | + |
| 28 | + gamma = rp.get_param("eos.gamma") |
| 29 | + |
| 30 | + # scale_height = rp.get_param("gresho.scale_height") |
| 31 | + dens_base = rp.get_param("gresho.dens_base") |
| 32 | + dens_cutoff = rp.get_param("gresho.dens_cutoff") |
| 33 | + |
| 34 | + rr = rp.get_param("gresho.r") |
| 35 | + u0 = rp.get_param("gresho.u0") |
| 36 | + p0 = rp.get_param("gresho.p0") |
| 37 | + |
| 38 | + # initialize the components -- we'll get a psure too |
| 39 | + # but that is used only to initialize the base state |
| 40 | + xmom[:, :] = 0.0 |
| 41 | + ymom[:, :] = 0.0 |
| 42 | + dens[:, :] = dens_cutoff |
| 43 | + |
| 44 | + # set the density to be stratified in the y-direction |
| 45 | + myg = my_data.grid |
| 46 | + pres = myg.scratch_array() |
| 47 | + |
| 48 | + dens[:, :] = dens_base |
| 49 | + |
| 50 | + pres[:, :] = p0 |
| 51 | + |
| 52 | + x_centre = 0.5 * (myg.x[0] + myg.x[-1]) |
| 53 | + y_centre = 0.5 * (myg.y[0] + myg.y[-1]) |
| 54 | + |
| 55 | + rad = numpy.sqrt((myg.x2d - x_centre)**2 + (myg.y2d - y_centre)**2) |
| 56 | + |
| 57 | + pres[rad <= rr] += 0.5 * (u0 * rad[rad <= rr]/rr)**2 |
| 58 | + pres[(rad > rr) & (rad <= 2*rr)] += \ |
| 59 | + u0**2 * (0.5 * (rad[(rad > rr) & (rad <= 2*rr)] / rr)**2 + |
| 60 | + 4 * (1 - rad[(rad > rr) & (rad <= 2*rr)]/rr + |
| 61 | + numpy.log(rad[(rad > rr) & (rad <= 2*rr)]/rr))) |
| 62 | + pres[rad > 2*rr] += u0**2 * (4 * numpy.log(2) - 2) |
| 63 | + # |
| 64 | + uphi = numpy.zeros_like(pres) |
| 65 | + uphi[rad <= rr] = u0 * rad[rad <= rr]/rr |
| 66 | + uphi[(rad > rr) & (rad <= 2*rr)] = u0 * (2 - rad[(rad > rr) & (rad <= 2*rr)]/rr) |
| 67 | + |
| 68 | + xmom[:, :] = -dens[:, :] * uphi[:, :] * (myg.y2d - y_centre) / rad[:, :] |
| 69 | + ymom[:, :] = dens[:, :] * uphi[:, :] * (myg.x2d - x_centre) / rad[:, :] |
| 70 | + |
| 71 | + ener[:, :] = pres[:, :]/(gamma - 1.0) + \ |
| 72 | + 0.5*(xmom[:, :]**2 + ymom[:, :]**2)/dens[:, :] |
| 73 | + |
| 74 | + eint = pres[:, :]/(gamma - 1.0) |
| 75 | + |
| 76 | + dens[:, :] = pres[:, :]/(eint[:, :]*(gamma - 1.0)) |
| 77 | + |
| 78 | + |
| 79 | +def finalize(): |
| 80 | + """ print out any information to the userad at the end of the run """ |
| 81 | + pass |
0 commit comments