diff --git a/compressible/interface.py b/compressible/interface.py index 7862ea807..bd0bdd194 100644 --- a/compressible/interface.py +++ b/compressible/interface.py @@ -419,7 +419,7 @@ def riemann_cgf(idir, ng, lambdastar_r = ustar + cstar_r if (pstar > p_r): - # the wave if a shock -- find the shock speed + # the wave is a shock -- find the shock speed sigma = (lambda_r + lambdastar_r) / 2.0 if (sigma > 0.0): @@ -713,7 +713,7 @@ def riemann_prim(idir, ng, lambdastar_r = ustar + cstar_r if (pstar > p_r): - # the wave if a shock -- find the shock speed + # the wave is a shock -- find the shock speed sigma = (lambda_r + lambdastar_r) / 2.0 if (sigma > 0.0): @@ -801,7 +801,7 @@ def riemann_prim(idir, ng, return q_int -@njit(cache=True) +# @njit(cache=True) def riemann_hllc(idir, ng, idens, ixmom, iymom, iener, irhoX, nspec, lower_solid, upper_solid, @@ -857,6 +857,9 @@ def riemann_hllc(idir, ng, # primitive variable states rho_l = U_l[i, j, idens] + if (rho_l <= 0): + print(f"-ve density {rho_l} at i,j = ({i},{j})") + # un = normal velocity; ut = transverse velocity if (idir == 1): un_l = U_l[i, j, ixmom] / rho_l diff --git a/compressible/problems/acoustic_pulse.py b/compressible/problems/acoustic_pulse.py index 1f98e5ada..aa36a9780 100644 --- a/compressible/problems/acoustic_pulse.py +++ b/compressible/problems/acoustic_pulse.py @@ -41,18 +41,19 @@ def init_data(myd, rp): ymin = rp.get_param("mesh.ymin") ymax = rp.get_param("mesh.ymax") - xctr = 0.5*(xmin + xmax) - yctr = 0.5*(ymin + ymax) + xctr = 0.5 * (xmin + xmax) + yctr = 0.5 * (ymin + ymax) dist = np.sqrt((myd.grid.x2d - xctr)**2 + (myd.grid.y2d - yctr)**2) dens[:, :] = rho0 idx = dist <= 0.5 - dens[idx] = rho0 + drho0*np.exp(-16*dist[idx]**2) * np.cos(np.pi*dist[idx])**6 + dens[idx] = rho0 + drho0 * \ + np.exp(-16 * dist[idx]**2) * np.cos(np.pi * dist[idx])**6 - p = (dens/rho0)**gamma - ener[:, :] = p/(gamma - 1) + p = (dens / rho0)**gamma + ener[:, :] = p / (gamma - 1) def finalize(): diff --git a/compressible/problems/advect.py b/compressible/problems/advect.py index 491338c2f..625d2e29b 100644 --- a/compressible/problems/advect.py +++ b/compressible/problems/advect.py @@ -38,22 +38,23 @@ def init_data(my_data, rp): ymin = rp.get_param("mesh.ymin") ymax = rp.get_param("mesh.ymax") - xctr = 0.5*(xmin + xmax) - yctr = 0.5*(ymin + ymax) + xctr = 0.5 * (xmin + xmax) + yctr = 0.5 * (ymin + ymax) # this is identical to the advection/smooth problem - dens[:, :] = 1.0 + np.exp(-60.0*((my_data.grid.x2d-xctr)**2 + - (my_data.grid.y2d-yctr)**2)) + dens[:, :] = 1.0 + np.exp(-60.0 * ((my_data.grid.x2d - xctr)**2 + + (my_data.grid.y2d - yctr)**2)) # velocity is diagonal u = 1.0 v = 1.0 - xmom[:, :] = dens[:, :]*u - ymom[:, :] = dens[:, :]*v + xmom[:, :] = dens[:, :] * u + ymom[:, :] = dens[:, :] * v # pressure is constant p = 1.0 - ener[:, :] = p/(gamma - 1.0) + 0.5*(xmom[:, :]**2 + ymom[:, :]**2)/dens[:, :] + ener[:, :] = p / (gamma - 1.0) + 0.5 * (xmom[:, :] + ** 2 + ymom[:, :]**2) / dens[:, :] def finalize(): diff --git a/compressible/problems/test.py b/compressible/problems/test.py index 2e463f6ce..f0a7b757f 100644 --- a/compressible/problems/test.py +++ b/compressible/problems/test.py @@ -9,7 +9,7 @@ def init_data(my_data, rp): # make sure that we are passed a valid patch object if not isinstance(my_data, patch.CellCenterData2d): - print("ERROR: patch invalid in sedov.py") + print("ERROR: patch invalid in test.py") print(my_data.__class__) sys.exit() diff --git a/compressible/simulation.py b/compressible/simulation.py index 42169d293..dc20cfa49 100644 --- a/compressible/simulation.py +++ b/compressible/simulation.py @@ -10,18 +10,18 @@ import compressible.derives as derives import compressible.unsplit_fluxes as flx import mesh.boundary as bnd -from simulation_null import NullSimulation, grid_setup, bc_setup +from simulation_null import NullSimulation, NullVariables, grid_setup, bc_setup import util.plot_tools as plot_tools import particles.particles as particles -class Variables(object): +class Variables(NullVariables): """ a container class for easy access to the different compressible variable by an integer key """ - def __init__(self, myd): - self.nvar = len(myd.names) + + def initialize(self, myd): # conserved variables -- we set these when we initialize for # they match the CellCenterData2d object @@ -58,19 +58,19 @@ def cons_to_prim(U, gamma, ivars, myg): q = myg.scratch_array(nvar=ivars.nq) q[:, :, ivars.irho] = U[:, :, ivars.idens] - q[:, :, ivars.iu] = U[:, :, ivars.ixmom]/U[:, :, ivars.idens] - q[:, :, ivars.iv] = U[:, :, ivars.iymom]/U[:, :, ivars.idens] + q[:, :, ivars.iu] = U[:, :, ivars.ixmom] / U[:, :, ivars.idens] + q[:, :, ivars.iv] = U[:, :, ivars.iymom] / U[:, :, ivars.idens] e = (U[:, :, ivars.iener] - - 0.5*q[:, :, ivars.irho]*(q[:, :, ivars.iu]**2 + - q[:, :, ivars.iv]**2))/q[:, :, ivars.irho] + 0.5 * q[:, :, ivars.irho] * (q[:, :, ivars.iu]**2 + + q[:, :, ivars.iv]**2)) / q[:, :, ivars.irho] q[:, :, ivars.ip] = eos.pres(gamma, q[:, :, ivars.irho], e) if ivars.naux > 0: - for nq, nu in zip(range(ivars.ix, ivars.ix+ivars.naux), - range(ivars.irhox, ivars.irhox+ivars.naux)): - q[:, :, nq] = U[:, :, nu]/q[:, :, ivars.irho] + for nq, nu in zip(range(ivars.ix, ivars.ix + ivars.naux), + range(ivars.irhox, ivars.irhox + ivars.naux)): + q[:, :, nq] = U[:, :, nu] / q[:, :, ivars.irho] return q @@ -81,18 +81,18 @@ def prim_to_cons(q, gamma, ivars, myg): U = myg.scratch_array(nvar=ivars.nvar) U[:, :, ivars.idens] = q[:, :, ivars.irho] - U[:, :, ivars.ixmom] = q[:, :, ivars.iu]*U[:, :, ivars.idens] - U[:, :, ivars.iymom] = q[:, :, ivars.iv]*U[:, :, ivars.idens] + U[:, :, ivars.ixmom] = q[:, :, ivars.iu] * U[:, :, ivars.idens] + U[:, :, ivars.iymom] = q[:, :, ivars.iv] * U[:, :, ivars.idens] rhoe = eos.rhoe(gamma, q[:, :, ivars.ip]) - U[:, :, ivars.iener] = rhoe + 0.5*q[:, :, ivars.irho]*(q[:, :, ivars.iu]**2 + - q[:, :, ivars.iv]**2) + U[:, :, ivars.iener] = rhoe + 0.5 * q[:, :, ivars.irho] * (q[:, :, ivars.iu]**2 + + q[:, :, ivars.iv]**2) if ivars.naux > 0: - for nq, nu in zip(range(ivars.ix, ivars.ix+ivars.naux), - range(ivars.irhox, ivars.irhox+ivars.naux)): - U[:, :, nu] = q[:, :, nq]*q[:, :, ivars.irho] + for nq, nu in zip(range(ivars.ix, ivars.ix + ivars.naux), + range(ivars.irhox, ivars.irhox + ivars.naux)): + U[:, :, nu] = q[:, :, nq] * q[:, :, ivars.irho] return U @@ -113,7 +113,8 @@ def initialize(self, extra_vars=None, ng=4): # define solver specific boundary condition routines bnd.define_bc("hse", BC.user, is_solid=False) - bnd.define_bc("ramp", BC.user, is_solid=False) # for double mach reflection problem + # for double mach reflection problem + bnd.define_bc("ramp", BC.user, is_solid=False) bc, bc_xodd, bc_yodd = bc_setup(self.rp) @@ -182,10 +183,10 @@ def method_compute_timestep(self): u, v, cs = self.cc_data.get_var(["velocity", "soundspeed"]) # the timestep is min(dx/(|u| + cs), dy/(|v| + cs)) - xtmp = self.cc_data.grid.dx/(abs(u) + cs) - ytmp = self.cc_data.grid.dy/(abs(v) + cs) + xtmp = self.cc_data.grid.dx / (abs(u) + cs) + ytmp = self.cc_data.grid.dy / (abs(v) + cs) - self.dt = cfl*float(min(xtmp.min(), ytmp.min())) + self.dt = cfl * float(min(xtmp.min(), ytmp.min())) def evolve(self): """ @@ -211,19 +212,19 @@ def evolve(self): old_ymom = ymom.copy() # conservative update - dtdx = self.dt/myg.dx - dtdy = self.dt/myg.dy + dtdx = self.dt / myg.dx + dtdy = self.dt / myg.dy for n in range(self.ivars.nvar): var = self.cc_data.get_var_by_index(n) var.v()[:, :] += \ - dtdx*(Flux_x.v(n=n) - Flux_x.ip(1, n=n)) + \ - dtdy*(Flux_y.v(n=n) - Flux_y.jp(1, n=n)) + dtdx * (Flux_x.v(n=n) - Flux_x.ip(1, n=n)) + \ + dtdy * (Flux_y.v(n=n) - Flux_y.jp(1, n=n)) # gravitational source terms - ymom[:, :] += 0.5*self.dt*(dens[:, :] + old_dens[:, :])*grav - ener[:, :] += 0.5*self.dt*(ymom[:, :] + old_ymom[:, :])*grav + ymom[:, :] += 0.5 * self.dt * (dens[:, :] + old_dens[:, :]) * grav + ener[:, :] += 0.5 * self.dt * (ymom[:, :] + old_ymom[:, :]) * grav if self.particles is not None: self.particles.update_particles(self.dt) @@ -257,7 +258,7 @@ def dovis(self): u = q[:, :, ivars.iu] v = q[:, :, ivars.iv] p = q[:, :, ivars.ip] - e = eos.rhoe(gamma, p)/rho + e = eos.rhoe(gamma, p) / rho magvel = np.sqrt(u**2 + v**2) @@ -297,7 +298,7 @@ def dovis(self): # plot particles ax.scatter(particle_positions[:, 0], - particle_positions[:, 1], s=5, c=colors, alpha=0.8, cmap="Greys") + particle_positions[:, 1], s=5, c=colors, alpha=0.8, cmap="Greys") ax.set_xlim([myg.xmin, myg.xmax]) ax.set_ylim([myg.ymin, myg.ymax]) diff --git a/compressible_mapped/__init__.py b/compressible_mapped/__init__.py new file mode 100644 index 000000000..8e9910c89 --- /dev/null +++ b/compressible_mapped/__init__.py @@ -0,0 +1,8 @@ +"""A method-of-lines compressible hydrodynamics solver. Piecewise +constant reconstruction is done in space and a Runge-Kutta time +integration is used to advance the solutiion. + +""" +__all__ = ["simulation"] + +from .simulation import Simulation diff --git a/compressible_mapped/_defaults b/compressible_mapped/_defaults new file mode 100644 index 000000000..526658669 --- /dev/null +++ b/compressible_mapped/_defaults @@ -0,0 +1,24 @@ +[driver] +cfl = 0.8 + + +[eos] +gamma = 1.4 ; pres = rho ener (gamma - 1) + + +[compressible] +use_flattening = 1 ; apply flattening at shocks (1) + +z0 = 0.75 ; flattening z0 parameter +z1 = 0.85 ; flattening z1 parameter +delta = 0.33 ; flattening delta parameter + +cvisc = 0.1 ; artifical viscosity coefficient + +limiter = 2 ; limiter (0 = none, 1 = 2nd order, 2 = 4th order) + +temporal_method = RK4 ; integration method (see mesh/integration.py) + +grav = 0.0 ; gravitational acceleration (in y-direction) + +riemann = HLLC ; HLLC or CGF diff --git a/compressible_mapped/fluxes.py b/compressible_mapped/fluxes.py new file mode 100644 index 000000000..d3646c73f --- /dev/null +++ b/compressible_mapped/fluxes.py @@ -0,0 +1,169 @@ +""" +This is a 2nd-order PLM method for a method-of-lines integration +(i.e., no characteristic tracing). + +We wish to solve + +.. math:: + + U_t + F^x_x + F^y_y = H + +we want U_{i+1/2} -- the interface values that are input to +the Riemann problem through the faces for each zone. + +Taylor expanding *in space only* yields:: + + dU + U = U + 0.5 dx -- + i+1/2,j,L i,j dx + +""" + +import numpy as np + +import compressible.interface as interface +import compressible as comp +import mesh.reconstruction as reconstruction +import mesh.array_indexer as ai + +from util import msg + + +def fluxes(my_data, rp, ivars, solid, tc): + """ + unsplitFluxes returns the fluxes through the x and y interfaces by + doing an unsplit reconstruction of the interface values and then + solving the Riemann problem through all the interfaces at once + + currently we assume a gamma-law EOS + + We use the method of Calhoun, D. A., Helzel, C., & LeVeque, R. J. (2008). SIAM review, 50(4), 723-752 for the mapped grids. + + Parameters + ---------- + my_data : CellCenterData2d object + The data object containing the grid and advective scalar that + we are advecting. + rp : RuntimeParameters object + The runtime parameters for the simulation + vars : Variables object + The Variables object that tells us which indices refer to which + variables + tc : TimerCollection object + The timers we are using to profile + + Returns + ------- + out : ndarray, ndarray + The fluxes on the x- and y-interfaces + + """ + + tm_flux = tc.timer("unsplitFluxes") + tm_flux.begin() + + myg = my_data.grid + U = my_data.data + gamma = rp.get_param("eos.gamma") + + # ========================================================================= + # transform by rotating the velocity components + # ========================================================================= + + U_xl = myg.scratch_array(nvar=ivars.nvar) + U_xr = myg.scratch_array(nvar=ivars.nvar) + + U_yl = myg.scratch_array(nvar=ivars.nvar) + U_yr = myg.scratch_array(nvar=ivars.nvar) + + b = 3 + + nx, ny, _ = np.shape(U_xl.v(buf=b)) + + for i in range(nx): + for j in range(ny): + U_xl.v(buf=b)[i, j] = (my_data.R_fcx.v(buf=b)[i, j] @ + U.ip(-1, buf=b)[i, j]) + U_xr.v(buf=b)[i, j] = (my_data.R_fcx.v(buf=b)[i, j] @ + U.v(buf=b)[i, j]) + U_yl.v(buf=b)[i, j] = (my_data.R_fcy.v(buf=b)[i, j] @ + U.jp(-1, buf=b)[i, j]) + U_yr.v(buf=b)[i, j] = (my_data.R_fcy.v(buf=b)[i, j] @ + U.v(buf=b)[i, j]) + + # ========================================================================= + # construct the fluxes normal to the interfaces + # ========================================================================= + tm_riem = tc.timer("Riemann") + tm_riem.begin() + + riemann = rp.get_param("compressible.riemann") + + if riemann == "HLLC": + riemannFunc = interface.riemann_hllc + elif riemann == "CGF": + riemannFunc = interface.riemann_cgf + else: + msg.fail("ERROR: Riemann solver undefined") + + _fx = riemannFunc(1, myg.ng, + ivars.idens, ivars.ixmom, ivars.iymom, ivars.iener, ivars.irhox, ivars.naux, + solid.xl, solid.xr, + gamma, U_xl, U_xr) + + _fy = riemannFunc(2, myg.ng, + ivars.idens, ivars.ixmom, ivars.iymom, ivars.iener, ivars.irhox, ivars.naux, + solid.yl, solid.yr, + gamma, U_yl, U_yr) + + F_x = ai.ArrayIndexer(d=_fx, grid=myg) + F_y = ai.ArrayIndexer(d=_fy, grid=myg) + + tm_riem.end() + + # ========================================================================= + # transform back + # ========================================================================= + + nx = myg.qx - 2 * myg.ng + ny = myg.qy - 2 * myg.ng + ilo = myg.ng + ihi = myg.ng + nx + jlo = myg.ng + jhi = myg.ng + ny + + A_plus_delta_q_x = myg.scratch_array(nvar=ivars.nvar) + A_minus_delta_q_x = myg.scratch_array(nvar=ivars.nvar) + A_plus_delta_q_y = myg.scratch_array(nvar=ivars.nvar) + A_minus_delta_q_y = myg.scratch_array(nvar=ivars.nvar) + + for i in range(ilo - 1, ihi + 1): + for j in range(jlo - 1, jhi + 1): + + A_plus_delta_q_x[i, j, :] = myg.gamma_fcx[i, j] * \ + (my_data.R_fcx[i, j]).T @ \ + (interface.consFlux(1, gamma, ivars.idens, ivars.ixmom, ivars.iymom, + ivars.iener, ivars.irhox, ivars.naux, + U_xr[i, j, :]) - F_x[i, j, :]) + A_minus_delta_q_x[i, j, :] = myg.gamma_fcx[i, j] * \ + (my_data.R_fcx[i, j]).T @ \ + (F_x[i, j, :] - interface.consFlux(1, gamma, ivars.idens, ivars.ixmom, + ivars.iymom, ivars.iener, + ivars.irhox, ivars.naux, + U_xl[i, j, :])) + + A_plus_delta_q_y[i, j, :] = myg.gamma_fcy[i, j] * \ + (my_data.R_fcy[i, j]).T @ \ + (interface.consFlux(2, gamma, ivars.idens, ivars.ixmom, ivars.iymom, + ivars.iener, ivars.irhox, ivars.naux, + U_yr[i, j, :]) - F_y[i, j, :]) + A_minus_delta_q_y[i, j, :] = myg.gamma_fcy[i, j] * \ + (my_data.R_fcy[i, j]).T @ \ + (F_y[i, j, :] - interface.consFlux(2, gamma, ivars.idens, ivars.ixmom, + ivars.iymom, ivars.iener, + ivars.irhox, ivars.naux, + U_yl[i, j, :])) + + tm_flux.end() + + return A_plus_delta_q_x, A_minus_delta_q_x, A_plus_delta_q_y, A_minus_delta_q_y diff --git a/compressible_mapped/problems/__init__.py b/compressible_mapped/problems/__init__.py new file mode 100644 index 000000000..fd7367b35 --- /dev/null +++ b/compressible_mapped/problems/__init__.py @@ -0,0 +1,2 @@ +__all__ = ['acoustic_pulse', 'advect', 'bubble', 'hse', 'kh', 'logo', 'quad', 'rt', 'rt2', 'sedov', 'sod', + 'test', 'ramp'] diff --git a/compressible_mapped/problems/_acoustic_pulse.defaults b/compressible_mapped/problems/_acoustic_pulse.defaults new file mode 100644 index 000000000..e1f6fd0ff --- /dev/null +++ b/compressible_mapped/problems/_acoustic_pulse.defaults @@ -0,0 +1,5 @@ +[acoustic_pulse] +rho0 = 1.4 +drho0 = 0.14 + + diff --git a/compressible_mapped/problems/_advect.defaults b/compressible_mapped/problems/_advect.defaults new file mode 100644 index 000000000..726cc5844 --- /dev/null +++ b/compressible_mapped/problems/_advect.defaults @@ -0,0 +1,3 @@ +[advect] + + diff --git a/compressible_mapped/problems/_kh.defaults b/compressible_mapped/problems/_kh.defaults new file mode 100644 index 000000000..d2709ebe0 --- /dev/null +++ b/compressible_mapped/problems/_kh.defaults @@ -0,0 +1,8 @@ +[kh] +rho_1 = 1 +v_1 = -1.0 + +rho_2 = 2 +v_2 = 1.0 + + diff --git a/compressible_mapped/problems/_sod.defaults b/compressible_mapped/problems/_sod.defaults new file mode 100644 index 000000000..80b86ad80 --- /dev/null +++ b/compressible_mapped/problems/_sod.defaults @@ -0,0 +1,13 @@ +[sod] + +direction = x ; direction of the flow + +dens_left = 1.0 +dens_right = 0.125 + +u_left = 0.0 +u_right = 0.0 + +p_left = 1.0 +p_right = 0.1 + diff --git a/compressible_mapped/problems/acoustic_pulse.py b/compressible_mapped/problems/acoustic_pulse.py new file mode 100644 index 000000000..e1ac4d95d --- /dev/null +++ b/compressible_mapped/problems/acoustic_pulse.py @@ -0,0 +1,76 @@ +from __future__ import print_function + +import sys +import numpy as np +import sympy +from sympy.abc import x, y + +from util import msg +import mesh.mapped as mapped + + +def init_data(myd, rp): + """initialize the acoustic_pulse problem. This comes from + McCourquodale & Coella 2011""" + + msg.bold("initializing the acoustic pulse problem...") + + # make sure that we are passed a valid patch object + if not isinstance(myd, mapped.MappedCellCenterData2d): + print("ERROR: patch invalid in acoustic_pulse.py") + print(myd.__class__) + sys.exit() + + # get the density, momenta, and energy as separate variables + dens = myd.get_var("density") + xmom = myd.get_var("x-momentum") + ymom = myd.get_var("y-momentum") + ener = myd.get_var("energy") + + # initialize the components, remember, that ener here is rho*eint + # + 0.5*rho*v**2, where eint is the specific internal energy + # (erg/g) + xmom[:, :] = 0.0 + ymom[:, :] = 0.0 + + gamma = rp.get_param("eos.gamma") + + rho0 = rp.get_param("acoustic_pulse.rho0") + drho0 = rp.get_param("acoustic_pulse.drho0") + + xmin = rp.get_param("mesh.xmin") + xmax = rp.get_param("mesh.xmax") + + ymin = rp.get_param("mesh.ymin") + ymax = rp.get_param("mesh.ymax") + + myg = myd.grid + + X, Y = myg.physical_coords() + + # rather than map the center of the coordinate grid, find the center by + # locating the center of a diagonal across the physical domain + xctr_t, yctr_t = 0.5 * \ + (myg.physical_coords(xmin, ymin) + myg.physical_coords(xmax, ymax)) + + dist = np.sqrt((X - xctr_t)**2 + (Y - yctr_t)**2) + + dens[:, :] = rho0 + idx = dist <= 0.5 + dens[idx] = rho0 + drho0 * \ + np.exp(-16 * dist[idx]**2) * np.cos(np.pi * dist[idx])**6 + + p = (dens / rho0)**gamma + ener[:, :] = p / (gamma - 1) + + +def sym_map(myg): + + return sympy.Matrix([y**2, -2 * x**2 - y]) + + # return sympy.Matrix([x * sympy.cos(y), x * sympy.sin(y)]) + + +def finalize(): + """ print out any information to the user at the end of the run """ + pass diff --git a/compressible_mapped/problems/advect.py b/compressible_mapped/problems/advect.py new file mode 100644 index 000000000..e77126544 --- /dev/null +++ b/compressible_mapped/problems/advect.py @@ -0,0 +1,78 @@ +from __future__ import print_function + +import sys +import numpy as np +import sympy +from sympy.abc import x, y + +import mesh.mapped as mapped +from util import msg + + +def init_data(my_data, rp): + """ initialize a smooth advection problem for testing convergence """ + + msg.bold("initializing the advect problem...") + + # make sure that we are passed a valid patch object + if not isinstance(my_data, mapped.MappedCellCenterData2d): + print("ERROR: patch invalid in advect.py") + print(my_data.__class__) + sys.exit() + + # get the density, momenta, and energy as separate variables + dens = my_data.get_var("density") + xmom = my_data.get_var("x-momentum") + ymom = my_data.get_var("y-momentum") + ener = my_data.get_var("energy") + + # initialize the components, remember, that ener here is rho*eint + # + 0.5*rho*v**2, where eint is the specific internal energy + # (erg/g) + dens[:, :] = 1.0 + xmom[:, :] = 0.0 + ymom[:, :] = 0.0 + + gamma = rp.get_param("eos.gamma") + + xmin = rp.get_param("mesh.xmin") + xmax = rp.get_param("mesh.xmax") + + ymin = rp.get_param("mesh.ymin") + ymax = rp.get_param("mesh.ymax") + + myg = my_data.grid + + X, Y = myg.physical_coords() + + xctr, yctr = 0.5 * \ + (myg.physical_coords(xmin, ymin) + myg.physical_coords(xmax, ymax)) + + # this is identical to the advection/smooth problem + dens[:, :] = 1.0 + np.exp(-60.0 * ((X - xctr)**2 + + (Y - yctr)**2)) + + # velocity is diagonal + u = 1.0 + v = 1.0 + xmom[:, :] = dens[:, :] * u + ymom[:, :] = dens[:, :] * v + + # pressure is constant + p = 1.0 + ener[:, :] = p / (gamma - 1.0) + 0.5 * (xmom[:, :] + ** 2 + ymom[:, :]**2) / dens[:, :] + + +def sym_map(myg): + + return sympy.Matrix([0.5*x + y, -x + 3]) + + +def finalize(): + """ print out any information to the user at the end of the run """ + + msg = """ + """ + + print(msg) diff --git a/compressible_mapped/problems/inputs.acoustic_pulse b/compressible_mapped/problems/inputs.acoustic_pulse new file mode 100644 index 000000000..21a1f9485 --- /dev/null +++ b/compressible_mapped/problems/inputs.acoustic_pulse @@ -0,0 +1,37 @@ +[driver] +max_steps = 5000 +tmax = 0.24 +# 0.192 * dx +fix_dt = 1.5e-3 + +[compressible] +cvisc = 0.1 + + +[io] +basename = acoustic_pulse_ +dt_out = 0.03 + +[eos] +gamma = 1.4 + + +[mesh] +nx = 128 +ny = 128 +xmax = 1.0 +ymax = 1.0 + +xlboundary = outflow +xrboundary = outflow + +ylboundary = outflow +yrboundary = outflow + + +[acoustic_pulse] +rho0 = 1.4 +drho0 = 0.14 + +[vis] +dovis = 1 diff --git a/compressible_mapped/problems/inputs.acoustic_pulse.256 b/compressible_mapped/problems/inputs.acoustic_pulse.256 new file mode 100644 index 000000000..e0dc73e8a --- /dev/null +++ b/compressible_mapped/problems/inputs.acoustic_pulse.256 @@ -0,0 +1,38 @@ +[driver] +max_steps = 5000 +tmax = 0.24 +# 0.192 * dx +fix_dt = 7.5e-4 + + +[compressible] +cvisc = 0.1 + + +[io] +basename = acoustic_pulse_256_ +dt_out = 0.03 + + +[eos] +gamma = 1.4 + + +[mesh] +nx = 256 +ny = 256 +xmax = 1.0 +ymax = 1.0 + +xlboundary = outflow +xrboundary = outflow + +ylboundary = outflow +yrboundary = outflow + +[acoustic_pulse] +rho0 = 1.4 +drho0 = 0.14 + +[vis] +dovis = 1 diff --git a/compressible_mapped/problems/inputs.acoustic_pulse.512 b/compressible_mapped/problems/inputs.acoustic_pulse.512 new file mode 100644 index 000000000..3aae69719 --- /dev/null +++ b/compressible_mapped/problems/inputs.acoustic_pulse.512 @@ -0,0 +1,39 @@ +[driver] +max_steps = 5000 +tmax = 0.24 +# 0.192 * dx +fix_dt = 3.75e-4 + + +[compressible] +cvisc = 0.1 + + +[io] +basename = acoustic_pulse_512_ +dt_out = 0.03 + + +[eos] +gamma = 1.4 + + +[mesh] +nx = 512 +ny = 512 +xmax = 1.0 +ymax = 1.0 + +xlboundary = outflow +xrboundary = outflow + +ylboundary = outflow +yrboundary = outflow + + +[acoustic_pulse] +rho0 = 1.4 +drho0 = 0.14 + +[vis] +dovis = 1 diff --git a/compressible_mapped/problems/inputs.acoustic_pulse.64 b/compressible_mapped/problems/inputs.acoustic_pulse.64 new file mode 100644 index 000000000..295af6d17 --- /dev/null +++ b/compressible_mapped/problems/inputs.acoustic_pulse.64 @@ -0,0 +1,38 @@ +[driver] +max_steps = 5000 +tmax = 0.24 +# 0.192 * dx +fix_dt = 3.0e-3 + +[compressible] +cvisc = 0.1 +limiter = 0 + + +[io] +basename = acoustic_pulse_64_ +dt_out = 0.03 + +[eos] +gamma = 1.4 + + +[mesh] +nx = 64 +ny = 64 +xmax = 1.0 +ymax = 1.0 + +xlboundary = outflow +xrboundary = outflow + +ylboundary = outflow +yrboundary = outflow + + +[acoustic_pulse] +rho0 = 1.4 +drho0 = 0.14 + +[vis] +dovis = 1 diff --git a/compressible_mapped/problems/inputs.advect.128 b/compressible_mapped/problems/inputs.advect.128 new file mode 100644 index 000000000..f0def69b5 --- /dev/null +++ b/compressible_mapped/problems/inputs.advect.128 @@ -0,0 +1,36 @@ +[driver] +max_steps = 500 +tmax = 1.0 + +init_tstep_factor = 1.0 +fix_dt = 0.0025 + + +[compressible] +limiter = 0 +cvisc = 0.1 + + +[io] +basename = advect_ + + +[eos] +gamma = 1.4 + + +[mesh] +nx = 128 +ny = 128 +xmax = 1.0 +ymax = 1.0 + +xlboundary = outflow +xrboundary = outflow + +ylboundary = outflow +yrboundary = outflow + + +[vis] +dovis = 1 diff --git a/compressible_mapped/problems/inputs.advect.256 b/compressible_mapped/problems/inputs.advect.256 new file mode 100644 index 000000000..30b083840 --- /dev/null +++ b/compressible_mapped/problems/inputs.advect.256 @@ -0,0 +1,36 @@ +[driver] +max_steps = 1000 +tmax = 1.0 + +init_tstep_factor = 1.0 +fix_dt = 0.00125 + + +[compressible] +limiter = 0 +cvisc = 0.1 + + +[io] +basename = advect_ + + +[eos] +gamma = 1.4 + + +[mesh] +nx = 256 +ny = 256 +xmax = 1.0 +ymax = 1.0 + +xlboundary = outflow +xrboundary = outflow + +ylboundary = outflow +yrboundary = outflow + + +[vis] +dovis = 1 diff --git a/compressible_mapped/problems/inputs.advect.64 b/compressible_mapped/problems/inputs.advect.64 new file mode 100644 index 000000000..30cfdc59d --- /dev/null +++ b/compressible_mapped/problems/inputs.advect.64 @@ -0,0 +1,36 @@ +[driver] +max_steps = 500 +tmax = 1.0 + +init_tstep_factor = 1.0 +fix_dt = 0.005 + + +[compressible] +limiter = 0 +cvisc = 0.1 + + +[io] +basename = advect_ + + +[eos] +gamma = 1.4 + + +[mesh] +nx = 64 +ny = 64 +xmax = 1.0 +ymax = 1.0 + +xlboundary = outflow +xrboundary = outflow + +ylboundary = outflow +yrboundary = outflow + + +[vis] +dovis = 1 diff --git a/compressible_mapped/problems/inputs.kh b/compressible_mapped/problems/inputs.kh new file mode 100644 index 000000000..386094c28 --- /dev/null +++ b/compressible_mapped/problems/inputs.kh @@ -0,0 +1,43 @@ +[driver] +max_steps = 5000 +tmax = 2.0 +splitting = unsplit + + +[compressible] +limiter = 2 +cvisc = 0.1 + + +[io] +basename = kh_ +tplot = 0.01 + + +[eos] +gamma = 1.4 + + +[mesh] +nx = 64 +ny = 64 +xmax = 1.0 +ymax = 1.0 + +xlboundary = periodic +xrboundary = periodic + +ylboundary = periodic +yrboundary = periodic + + +[kh] +rho_1 = 1 +v_1 = -0.5 + +rho_2 = 2 +v_2 = 0.5 + + +[vis] +dovis = 1 diff --git a/compressible_mapped/problems/inputs.sod.x b/compressible_mapped/problems/inputs.sod.x new file mode 100644 index 000000000..bae148fdc --- /dev/null +++ b/compressible_mapped/problems/inputs.sod.x @@ -0,0 +1,32 @@ +# simple inputs files for the unsplit CTU hydro scheme + +[driver] +max_steps = 200 +tmax = 0.2 + +[compressible] +limiter = 1 + +[io] +basename = sod_x_ +dt_out = 0.05 + +[mesh] +nx = 128 +ny = 10 +xmax = 1.0 +ymax = .05 +xlboundary = outflow +xrboundary = outflow + +[sod] +direction = x + +dens_left = 1.0 +dens_right = 0.125 + +u_left = 0.0 +u_right = 0.0 + +p_left = 1.0 +p_right = 0.1 diff --git a/compressible_mapped/problems/inputs.sod.y b/compressible_mapped/problems/inputs.sod.y new file mode 100644 index 000000000..9e952fd5d --- /dev/null +++ b/compressible_mapped/problems/inputs.sod.y @@ -0,0 +1,29 @@ +# simple inputs files for the unsplit CTU hydro scheme + +[driver] +max_steps = 200 +tmax = 0.2 + +[io] +basename = sod_y_ +dt_out = 0.05 + +[mesh] +nx = 10 +ny = 128 +xmax = .05 +ymax = 1.0 +ylboundary = outflow +yrboundary = outflow + +[sod] +direction = y + +dens_left = 1.0 +dens_right = 0.125 + +u_left = 0.0 +u_right = 0.0 + +p_left = 1.0 +p_right = 0.1 diff --git a/compressible_mapped/problems/kh.py b/compressible_mapped/problems/kh.py new file mode 100644 index 000000000..6b8662a3b --- /dev/null +++ b/compressible_mapped/problems/kh.py @@ -0,0 +1,97 @@ +from __future__ import print_function + +import numpy as np +import sympy +from sympy.abc import x, y + +import mesh.mapped as mapped +from util import msg + + +def init_data(my_data, rp): + """ initialize the Kelvin-Helmholtz problem """ + + msg.bold("initializing the Kelvin-Helmholtz problem...") + + # make sure that we are passed a valid patch object + if not isinstance(my_data, mapped.MappedCellCenterData2d): + print(my_data.__class__) + msg.fail("ERROR: patch invalid in kh.py") + + # get the density, momenta, and energy as separate variables + dens = my_data.get_var("density") + xmom = my_data.get_var("x-momentum") + ymom = my_data.get_var("y-momentum") + ener = my_data.get_var("energy") + + # initialize the components, remember, that ener here is rho*eint + # + 0.5*rho*v**2, where eint is the specific internal energy + # (erg/g) + dens[:, :] = 1.0 + xmom[:, :] = 0.0 + ymom[:, :] = 0.0 + + rho_1 = rp.get_param("kh.rho_1") + v_1 = rp.get_param("kh.v_1") + rho_2 = rp.get_param("kh.rho_2") + v_2 = rp.get_param("kh.v_2") + + xmin_coord = rp.get_param("mesh.xmin") + ymin_coord = rp.get_param("mesh.ymin") + xmax_coord = rp.get_param("mesh.xmax") + ymax_coord = rp.get_param("mesh.ymax") + + gamma = rp.get_param("eos.gamma") + + myg = my_data.grid + + xmin, ymin = myg.physical_coords(xmin_coord, ymin_coord) + xmax, ymax = myg.physical_coords(xmax_coord, ymax_coord) + + dy = 0.025*(ymax - ymin) + w0 = 0.01 + vm = 0.5 * (v_1 - v_2) + rhom = 0.5 * (rho_1 - rho_2) + + X, Y = myg.physical_coords() + + idx1 = Y < 0.25*ymax + idx2 = np.logical_and(Y >= 0.25*ymax, Y < 0.5*ymax) + idx3 = np.logical_and(Y >= 0.5*ymax, Y < 0.75*ymax) + idx4 = Y >= 0.75*ymax + + # we will initialize momemum as velocity for now + + # lower quarter + dens[idx1] = rho_1 - rhom * np.exp((Y[idx1] - 0.25*ymax) / dy) + xmom[idx1] = v_1 - vm * np.exp((Y[idx1] - 0.25*ymax) / dy) + + # second quarter + dens[idx2] = rho_2 + rhom * np.exp((0.25*ymax - Y[idx2]) / dy) + xmom[idx2] = v_2 + vm * np.exp((0.25*ymax - Y[idx2]) / dy) + + # third quarter + dens[idx3] = rho_2 + rhom * np.exp((Y[idx3] - 0.75*ymax) / dy) + xmom[idx3] = v_2 + vm * np.exp((Y[idx3] - 0.75*ymax) / dy) + + # fourth quarter + dens[idx4] = rho_1 - rhom * np.exp((0.75*ymax - Y[idx4]) / dy) + xmom[idx4] = v_1 - vm * np.exp((0.75*ymax - Y[idx4]) / dy) + + # upper half + xmom[:, :] *= dens + ymom[:, :] = dens * w0 * np.sin(4 * np.pi * X/(xmax-xmin)) + + p = 2.5 + ener[:, :] = p / (gamma - 1.0) + 0.5 * (xmom[:, :] + ** 2 + ymom[:, :]**2) / dens[:, :] + + +def sym_map(myg): + + return sympy.Matrix([-x + 4, 2*y]) + + +def finalize(): + """ print out any information to the user at the end of the run """ + pass diff --git a/compressible_mapped/problems/sod.py b/compressible_mapped/problems/sod.py new file mode 100644 index 000000000..7c7b2241e --- /dev/null +++ b/compressible_mapped/problems/sod.py @@ -0,0 +1,109 @@ +from __future__ import print_function + +import sys +import sympy +from sympy.abc import x, y + +import mesh.mapped as mapped +from util import msg + + +def init_data(my_data, rp): + """ initialize the sod problem """ + + msg.bold("initializing the sod problem...") + + # make sure that we are passed a valid patch object + if not isinstance(my_data, mapped.MappedCellCenterData2d): + print("ERROR: patch invalid in sod.py") + print(my_data.__class__) + sys.exit() + + # get the sod parameters + dens_left = rp.get_param("sod.dens_left") + dens_right = rp.get_param("sod.dens_right") + + u_left = rp.get_param("sod.u_left") + u_right = rp.get_param("sod.u_right") + + p_left = rp.get_param("sod.p_left") + p_right = rp.get_param("sod.p_right") + + # get the density, momenta, and energy as separate variables + dens = my_data.get_var("density") + xmom = my_data.get_var("x-momentum") + ymom = my_data.get_var("y-momentum") + ener = my_data.get_var("energy") + + # initialize the components, remember, that ener here is rho*eint + # + 0.5*rho*v**2, where eint is the specific internal energy + # (erg/g) + xmin = rp.get_param("mesh.xmin") + xmax = rp.get_param("mesh.xmax") + + ymin = rp.get_param("mesh.ymin") + ymax = rp.get_param("mesh.ymax") + + gamma = rp.get_param("eos.gamma") + + direction = rp.get_param("sod.direction") + + myg = my_data.grid + + X, Y = myg.physical_coords() + + xctr, yctr = 0.5 * \ + (myg.physical_coords(xmin, ymin) + myg.physical_coords(xmax, ymax)) + + if direction == "x": + + # left + idxl = X <= xctr + + dens[idxl] = dens_left + xmom[idxl] = dens_left * u_left + ymom[idxl] = 0.0 + ener[idxl] = p_left / (gamma - 1.0) + 0.5 * xmom[idxl] * u_left + + # right + idxr = X > xctr + + dens[idxr] = dens_right + xmom[idxr] = dens_right * u_right + ymom[idxr] = 0.0 + ener[idxr] = p_right / (gamma - 1.0) + 0.5 * xmom[idxr] * u_right + + else: + + # bottom + idxb = Y <= yctr + + dens[idxb] = dens_left + xmom[idxb] = 0.0 + ymom[idxb] = dens_left * u_left + ener[idxb] = p_left / (gamma - 1.0) + 0.5 * ymom[idxb] * u_left + + # top + idxt = Y > yctr + + dens[idxt] = dens_right + xmom[idxt] = 0.0 + ymom[idxt] = dens_right * u_right + ener[idxt] = p_right / (gamma - 1.0) + 0.5 * ymom[idxt] * u_right + + +def sym_map(myg): + + return sympy.Matrix([x - y, x + y + 3]) + + +def finalize(): + """ print out any information to the user at the end of the run """ + + msg = """ + The script analysis/sod_compare.py can be used to compare + this output to the exact solution. Some sample exact solution + data is present as analysis/sod-exact.out + """ + + print(msg) diff --git a/compressible_mapped/problems/test.py b/compressible_mapped/problems/test.py new file mode 100644 index 000000000..035d74bf7 --- /dev/null +++ b/compressible_mapped/problems/test.py @@ -0,0 +1,41 @@ +from __future__ import print_function + +import sys +import sympy +from sympy.abc import x, y + +import mesh.mapped as mapped + + +def init_data(my_data, rp): + """ an init routine for unit testing """ + + # make sure that we are passed a valid patch object + if not isinstance(my_data, mapped.MappedCellCenterData2d): + print("ERROR: patch invalid in test.py") + print(my_data.__class__) + sys.exit() + + # get the density, momenta, and energy as separate variables + dens = my_data.get_var("density") + xmom = my_data.get_var("x-momentum") + ymom = my_data.get_var("y-momentum") + ener = my_data.get_var("energy") + + # initialize the components, remember, that ener here is rho*eint + # + 0.5*rho*v**2, where eint is the specific internal energy + # (erg/g) + dens[:, :] = 1.0 + xmom[:, :] = 0.0 + ymom[:, :] = 0.0 + ener[:, :] = 2.5 + + +def sym_map(myg): + + return sympy.Matrix([-2 * y, x + 4]) + + +def finalize(): + """ print out any information to the user at the end of the run """ + pass diff --git a/compressible_mapped/simulation.py b/compressible_mapped/simulation.py new file mode 100644 index 000000000..08e58f049 --- /dev/null +++ b/compressible_mapped/simulation.py @@ -0,0 +1,307 @@ +from __future__ import print_function + +import numpy as np +import importlib +import matplotlib.pyplot as plt + +import mesh.integration as integration +import compressible +import compressible_rk +import compressible_mapped.fluxes as flx +from util import msg, plot_tools +import mesh.mapped as mapped +import mesh.boundary as bnd +import compressible.BC as BC +import compressible.eos as eos +from simulation_null import bc_setup +import particles.particles as particles +import compressible.derives as derives + + +def mapped_grid_setup(rp, map, ng=1): + nx = rp.get_param("mesh.nx") + ny = rp.get_param("mesh.ny") + + try: + xmin = rp.get_param("mesh.xmin") + except KeyError: + xmin = 0.0 + msg.warning("mesh.xmin not set, defaulting to 0.0") + + try: + xmax = rp.get_param("mesh.xmax") + except KeyError: + xmax = 1.0 + msg.warning("mesh.xmax not set, defaulting to 1.0") + + try: + ymin = rp.get_param("mesh.ymin") + except KeyError: + ymin = 0.0 + msg.warning("mesh.ymin not set, defaulting to 0.0") + + try: + ymax = rp.get_param("mesh.ymax") + except KeyError: + ymax = 1.0 + msg.warning("mesh.ynax not set, defaulting to 1.0") + + my_grid = mapped.MappedGrid2d(map, nx, ny, + xmin=xmin, xmax=xmax, + ymin=ymin, ymax=ymax, ng=ng) + return my_grid + +def mapped_bc_setup(rp): + + bnd.define_bc("mapped_periodic", BC.user, is_solid=False) + bnd.define_bc("mapped_outflow", BC.user, is_solid=False) + bnd.define_bc("mapped_neumann", BC.user, is_solid=False) + bnd.define_bc("mapped_reflect", BC.user, is_solid=True) + bnd.define_bc("mapped_reflect-even", BC.user, is_solid=True) + bnd.define_bc("mapped_reflect-odd", BC.user, is_solid=True) + bnd.define_bc("mapped_dirichlet", BC.user, is_solid=True) + + # first figure out the BCs + try: + xlb_type = "mapped_" + rp.get_param("mesh.xlboundary") + except KeyError: + xlb_type = "mapped_periodic" + msg.warning("mesh.xlboundary is not set, defaulting to mapped_periodic") + + try: + xrb_type = "mapped_" + rp.get_param("mesh.xrboundary") + except KeyError: + xrb_type = "mapped_periodic" + msg.warning("mesh.xrboundary is not set, defaulting to mapped_periodic") + + try: + ylb_type = "mapped_" + rp.get_param("mesh.ylboundary") + except KeyError: + ylb_type = "mapped_periodic" + msg.warning("mesh.ylboundary is not set, defaulting to mapped_periodic") + + try: + yrb_type = "mapped_" + rp.get_param("mesh.yrboundary") + except KeyError: + yrb_type = "mapped_periodic" + msg.warning("mesh.yrboundary is not set, defaulting to mapped_periodic") + + bc = bnd.BC(xlb=xlb_type, xrb=xrb_type, + ylb=ylb_type, yrb=yrb_type) + + # if we are reflecting, we need odd reflection in the normal + # directions for the velocity + bc_xodd = bnd.BC(xlb=xlb_type, xrb=xrb_type, + ylb=ylb_type, yrb=yrb_type, + odd_reflect_dir="x") + + bc_yodd = bnd.BC(xlb=xlb_type, xrb=xrb_type, + ylb=ylb_type, yrb=yrb_type, + odd_reflect_dir="y") + + return bc, bc_xodd, bc_yodd + + +class Simulation(compressible_rk.Simulation): + """The main simulation class for the method of lines compressible + hydrodynamics solver""" + + def __init__(self, solver_name, problem_name, rp, timers=None): + + super().__init__(solver_name, problem_name, rp, timers=timers, + data_class=mapped.MappedCellCenterData2d) + + def initialize(self, extra_vars=None, ng=4): + """ + Initialize the grid and variables for compressible flow and set + the initial conditions for the chosen problem. + """ + problem = importlib.import_module("{}.problems.{}".format( + self.solver_name, self.problem_name)) + + my_grid = mapped_grid_setup(self.rp, problem.sym_map, ng=ng) + my_data = self.data_class(my_grid) + + # # define solver specific boundary condition routines + # bnd.define_bc("hse", BC.user, is_solid=False) + # # for double mach reflection problem + # bnd.define_bc("ramp", BC.user, is_solid=False) + + bc, bc_xodd, bc_yodd = bc_setup(self.rp) + + # are we dealing with solid boundaries? we'll use these for + # the Riemann solver + self.solid = bnd.bc_is_solid(bc) + + # density and energy + my_data.register_var("density", bc) + my_data.register_var("energy", bc) + my_data.register_var("x-momentum", bc_xodd) + my_data.register_var("y-momentum", bc_yodd) + + # any extras? + if extra_vars is not None: + for v in extra_vars: + my_data.register_var(v, bc) + + # store the EOS gamma as an auxillary quantity so we can have a + # self-contained object stored in output files to make plots. + # store grav because we'll need that in some BCs + my_data.set_aux("gamma", self.rp.get_param("eos.gamma")) + + my_data.create() + + self.cc_data = my_data + + if self.rp.get_param("particles.do_particles") == 1: + self.particles = particles.Particles(self.cc_data, bc, self.rp) + + # some auxillary data that we'll need to fill GC in, but isn't + # really part of the main solution + aux_data = self.data_class(my_grid) + aux_data.register_var("ymom_src", bc_yodd) + aux_data.register_var("E_src", bc) + aux_data.create() + self.aux_data = aux_data + + self.ivars = compressible.Variables(my_data) + + # make rotation matrices + my_data.make_rotation_matrices(self.ivars) + + # derived variables + self.cc_data.add_derived(derives.derive_primitives) + + # initial conditions for the problem + problem.init_data(self.cc_data, self.rp) + + if self.verbose > 0: + print(my_data) + + def substep(self, myd): + """ + take a single substep in the RK timestepping starting with the + conservative state defined as part of myd + """ + + myg = myd.grid + k = myg.scratch_array(nvar=self.ivars.nvar) + + flux_xp, flux_xm, flux_yp, flux_ym = flx.fluxes(myd, self.rp, + self.ivars, self.solid, self.tc) + b = 3 + for n in range(self.ivars.nvar): + k.v(n=n, buf=b)[:, :] = \ + (flux_xp.v(n=n, buf=b) + flux_xm.ip(1, n=n, buf=b)) / (myg.dx * myg.kappa.v(buf=b)) + \ + (flux_yp.v(n=n, buf=b) + flux_ym.jp(1, n=n, buf=b)) / \ + (myg.dy * myg.kappa.v(buf=b)) + + return -k + + def evolve(self): + """ + Evolve the equations of compressible hydrodynamics through a + timestep dt. + """ + + tm_evolve = self.tc.timer("evolve") + tm_evolve.begin() + + myd = self.cc_data + + method = self.rp.get_param("compressible.temporal_method") + + rk = integration.RKIntegrator(myd.t, self.dt, method=method) + rk.set_start(myd) + + for s in range(rk.nstages()): + ytmp = rk.get_stage_start( + s, clone_function=mapped.mapped_cell_center_data_clone) + ytmp.fill_BC_all() + k = self.substep(ytmp) + rk.store_increment(s, k) + + rk.compute_final_update() + + if self.particles is not None: + self.particles.update_particles(self.dt) + + # increment the time + myd.t += self.dt + self.n += 1 + + tm_evolve.end() + + def dovis(self): + """ + Do runtime visualization. + """ + + plt.clf() + + plt.rc("font", size=10) + + # we do this even though ivars is in self, so this works when + # we are plotting from a file + ivars = compressible.Variables(self.cc_data) + + # access gamma from the cc_data object so we can use dovis + # outside of a running simulation. + gamma = self.cc_data.get_aux("gamma") + + q = compressible.cons_to_prim( + self.cc_data.data, gamma, ivars, self.cc_data.grid) + + rho = q[:, :, ivars.irho] + u = q[:, :, ivars.iu] + v = q[:, :, ivars.iv] + p = q[:, :, ivars.ip] + e = eos.rhoe(gamma, p) / rho + + magvel = np.sqrt(u**2 + v**2) + + myg = self.cc_data.grid + + fields = [rho, magvel, p, e] + field_names = [r"$\rho$", r"U", "p", "e"] + + _, axes, cbar_title = plot_tools.setup_axes(myg, len(fields)) + + X, Y = myg.physical_coords() + X = X[myg.ng:-myg.ng, myg.ng:-myg.ng] + Y = Y[myg.ng:-myg.ng, myg.ng:-myg.ng] + + for n, ax in enumerate(axes): + v = fields[n] + + img = ax.pcolormesh(X, Y, v.v(), cmap=self.cm) + + ax.set_xlabel("x") + ax.set_ylabel("y") + + # needed for PDF rendering + cb = axes.cbar_axes[n].colorbar(img) + cb.solids.set_rasterized(True) + cb.solids.set_edgecolor("face") + + if cbar_title: + cb.ax.set_title(field_names[n]) + else: + ax.set_title(field_names[n]) + + if self.particles is not None: + ax = axes[0] + particle_positions = self.particles.get_positions() + # dye particles + colors = self.particles.get_init_positions()[:, 0] + + # plot particles + ax.scatter(particle_positions[:, 0], + particle_positions[:, 1], s=5, c=colors, alpha=0.8, cmap="Greys") + ax.set_xlim([myg.xmin, myg.xmax]) + ax.set_ylim([myg.ymin, myg.ymax]) + + plt.figtext(0.05, 0.0125, "t = {:10.5g}".format(self.cc_data.t)) + + plt.pause(0.001) + plt.draw() diff --git a/compressible_mapped/tests/test_compressible_mapped.py b/compressible_mapped/tests/test_compressible_mapped.py new file mode 100644 index 000000000..c589b2b98 --- /dev/null +++ b/compressible_mapped/tests/test_compressible_mapped.py @@ -0,0 +1,66 @@ +import numpy as np +from numpy.testing import assert_array_equal +import pytest +import sympy + +from util import runparams +import compressible_mapped.simulation as sn + + +class TestSimulation(object): + @classmethod + def setup_class(cls): + """ this is run once for each class before any tests """ + pass + + @classmethod + def teardown_class(cls): + """ this is run once for each class after all tests """ + pass + + def setup_method(self): + """ this is run before each test """ + self.rp = runparams.RuntimeParameters() + + self.rp.params["mesh.nx"] = 8 + self.rp.params["mesh.ny"] = 8 + self.rp.params["particles.do_particles"] = 0 + + self.rp.params["eos.gamma"] = 1.4 + + self.sim = sn.Simulation("compressible_mapped", "test", self.rp) + self.sim.initialize() + + def teardown_method(self): + """ this is run after each test """ + self.rp = None + self.sim = None + + def test_initialization_state(self): + """ + Test state initialized properly + """ + dens = self.sim.cc_data.get_var("density") + assert dens.min() == 1.0 and dens.max() == 1.0 + + ener = self.sim.cc_data.get_var("energy") + assert ener.min() == 2.5 and ener.max() == 2.5 + + def test_mapped_grid(self): + + myg = self.sim.cc_data.grid + + kappa = myg.scratch_array() + 2 + assert_array_equal(kappa, myg.kappa) + + gamma_fcx = np.ones_like(kappa) * 2 + assert_array_equal(gamma_fcx, myg.gamma_fcx) + + gamma_fcy = np.ones_like(kappa) + assert_array_equal(gamma_fcy, myg.gamma_fcy) + + R_x = sympy.Matrix([[0, -1], [1, 0]]) + R_y = sympy.Matrix([[0, 1], [-1, 0]]) + + assert_array_equal(myg.R_fcx(2, 0, 1)[3, 5], R_x) + assert_array_equal(myg.R_fcy(2, 0, 1)[6, 0], R_y) diff --git a/docs/source/compressible_mapped.problems.rst b/docs/source/compressible_mapped.problems.rst new file mode 100644 index 000000000..79f378c45 --- /dev/null +++ b/docs/source/compressible_mapped.problems.rst @@ -0,0 +1,58 @@ +compressible\_mapped\.problems package +====================================== + +.. automodule:: compressible_mapped.problems + :members: + :undoc-members: + :show-inheritance: + +Submodules +---------- + +compressible\_mapped\.problems\.acoustic\_pulse module +------------------------------------------------------ + +.. automodule:: compressible_mapped.problems.acoustic_pulse + :members: + :undoc-members: + :show-inheritance: + +compressible\_mapped\.problems\.advect module +--------------------------------------------- + +.. automodule:: compressible_mapped.problems.advect + :members: + :undoc-members: + :show-inheritance: + +compressible\_mapped\.problems\.kh module +----------------------------------------- + +.. automodule:: compressible_mapped.problems.kh + :members: + :undoc-members: + :show-inheritance: + +compressible\_mapped\.problems\.sedov module +-------------------------------------------- + +.. automodule:: compressible_mapped.problems.sedov + :members: + :undoc-members: + :show-inheritance: + +compressible\_mapped\.problems\.sod module +------------------------------------------ + +.. automodule:: compressible_mapped.problems.sod + :members: + :undoc-members: + :show-inheritance: + +compressible\_mapped\.problems\.test module +------------------------------------------- + +.. automodule:: compressible_mapped.problems.test + :members: + :undoc-members: + :show-inheritance: diff --git a/docs/source/compressible_mapped.rst b/docs/source/compressible_mapped.rst new file mode 100644 index 000000000..e7deb8225 --- /dev/null +++ b/docs/source/compressible_mapped.rst @@ -0,0 +1,33 @@ +compressible\_mapped package +========================= + +.. automodule:: compressible_mapped + :members: + :undoc-members: + :show-inheritance: + +Subpackages +----------- + +.. toctree:: + + compressible_mapped.problems + +Submodules +---------- + +compressible\_mapped\.fluxes module +-------------------------------- + +.. automodule:: compressible_mapped.fluxes + :members: + :undoc-members: + :show-inheritance: + +compressible\_mapped\.simulation module +------------------------------------ + +.. automodule:: compressible_mapped.simulation + :members: + :undoc-members: + :show-inheritance: diff --git a/docs/source/compressible_mapped_defaults.inc b/docs/source/compressible_mapped_defaults.inc new file mode 100644 index 000000000..5417b0298 --- /dev/null +++ b/docs/source/compressible_mapped_defaults.inc @@ -0,0 +1,37 @@ +* section: [compressible] + + +----------------------------------+----------------+----------------------------------------------------+ + | option | value | description | + +==================================+================+====================================================+ + | use_flattening | ``1`` | apply flattening at shocks (1) | + +----------------------------------+----------------+----------------------------------------------------+ + | z0 | ``0.75`` | flattening z0 parameter | + +----------------------------------+----------------+----------------------------------------------------+ + | z1 | ``0.85`` | flattening z1 parameter | + +----------------------------------+----------------+----------------------------------------------------+ + | delta | ``0.33`` | flattening delta parameter | + +----------------------------------+----------------+----------------------------------------------------+ + | cvisc | ``0.1`` | artifical viscosity coefficient | + +----------------------------------+----------------+----------------------------------------------------+ + | limiter | ``2`` | limiter (0 = none, 1 = 2nd order, 2 = 4th order) | + +----------------------------------+----------------+----------------------------------------------------+ + | temporal_method | ``RK4`` | integration method (see mesh/integration.py) | + +----------------------------------+----------------+----------------------------------------------------+ + | riemann | ``HLLC`` | riemann solver (HLLC or CGF) | + +----------------------------------+----------------+----------------------------------------------------+ + +* section: [driver] + + +----------------------------------+----------------+----------------------------------------------------+ + | option | value | description | + +==================================+================+====================================================+ + | cfl | ``0.8`` | | + +----------------------------------+----------------+----------------------------------------------------+ + +* section: [eos] + + +----------------------------------+----------------+----------------------------------------------------+ + | option | value | description | + +==================================+================+====================================================+ + | gamma | ``1.4`` | pres = rho ener (gamma - 1) | + +----------------------------------+----------------+----------------------------------------------------+ diff --git a/docs/source/mapped_basics.rst b/docs/source/mapped_basics.rst new file mode 100644 index 000000000..3eaa0baca --- /dev/null +++ b/docs/source/mapped_basics.rst @@ -0,0 +1,82 @@ +:mod:`mesh.mapped ` implementation and use +------------------------------------------------------- + +This module defines the structures required for mapped grids. It defines the classes +:func:`MappedGrid2d ` and +:func:`MappedCellCenterData2d ` that inherit from +the un-mapped versions, :func:`Grid2d ` and +:func:`CellCenterData2d `. + +We import the basic mapped grids functionality as: + +.. code-block:: python + + import mesh.mapped as mapped + +The classes in the mapped module that we interact with are: + +* :func:`mapped.MappedGrid2d `: this is the main mapped + grid class. It holds the coordinate grid (a :func:`Grid2d ` + object) and details about the mapping from the physical grid to the + coordinate grid. The mapping is defined using a sympy function: when the object + is instantiated, this analytic expression with then be used to calculate + the various numerical quantities on the actual grid. + +* :func:`mapped.MappedCellCenterData2d `: this + holds the cell-centered data on the coordinate grid. It is identical to its + parent class, :func:`CellCenterData2d `, except + it also defines the rotation matrices that define how the velocity + components are mapped between the physical and coordinate grids (these + rotation matrices cannot be fully defined on the grid itself as they require + information about the variables which the grid alone does not have). + +The procedure for setting up a mapped grid and the data that lives on it is as +follows: + +.. code-block:: python + + import sympy + from sympy.abs import x, y + import mesh.mapped as mapped + import mesh.boundary as bnd + import compressible + + def mapping(myg): + return sympy.Matrix([2 * x + 3, y**2]) + + myg = mapped.MappedGrid2d(mapping, 16, 32, xmax=1.0, ymax=2.0) + +This creates a 2-d mapped grid object ``myg`` with 16 zones in the x-direction +and 32 zones in the y-direction on the coordinate grid. It also specifies the +extent of the coordinate grid in the x- and y-directions. The mapping is defined +using a ``sympy.Matrix`` object as a function of the cartesian coordinates +``x`` and ``y``. The mapping itself is a function of the mapped grid object. + +.. note:: + + The mapping must be specified as a function of the ``sympy.Symbol`` objects + ``x`` and ``y`` only. + +.. code-block:: python + + mydata = mapped.MappedCellCenterData2d(myg) + + bc = bnd.BC(xlb="periodic", xrb="periodic", ylb="reflect-even", yrb="outflow") + + mydata.register_var("a", bc) + mydata.create() + + ivars = compressible.Variables(mydata) + + mydata.make_rotation_matrices(ivars) + +This creates the mapped cell-centered data object, ``mydata``, that lives on the +mapped grid we just defined. Just as we would for the non-mapped grid, we next +create a boundary condition object to specify the boundary conditions on each edge, +register variables and call the ``create()`` function to allocate storage for the +variables. For the mapped grid, there is one more step that we must do: define the +rotation matrices associated with the variables. To do this, we define a +``Variables`` object (here we use the one from the compressible module), then +pass this to the data object's :func:`make_rotation_matrices ` method. This will +take the rotation matrix function defined on the mapped grid object ``myg`` and +create the actual matrices associated with the data's variables. diff --git a/docs/source/mesh.rst b/docs/source/mesh.rst index 0b8fde256..0503b2dee 100644 --- a/docs/source/mesh.rst +++ b/docs/source/mesh.rst @@ -41,6 +41,14 @@ mesh\.integration module :undoc-members: :show-inheritance: +mesh\.mapped module +------------------------ + +.. automodule:: mesh.mapped + :members: + :undoc-members: + :show-inheritance: + mesh\.patch module ------------------ diff --git a/docs/source/mesh_basics.rst b/docs/source/mesh_basics.rst index 5ba86e6f8..cc6b1ce29 100644 --- a/docs/source/mesh_basics.rst +++ b/docs/source/mesh_basics.rst @@ -27,7 +27,7 @@ We import the basic mesh functionality as: import mesh.boundary as bnd import mesh.array_indexer as ai -There are several main objects in the patch class that we interact with: +There are several main classes in the patch module that we interact with: * :func:`patch.Grid2d `: this is the main grid object. It is basically a container that holds the number of zones @@ -96,6 +96,8 @@ Note that each variable needs to specify a BC—this allows us to do different actions for each variable (for example, some may do even reflection while others may do odd reflection). +.. include:: mapped_basics.rst + Jupyter notebook ---------------- diff --git a/mesh/__init__.py b/mesh/__init__.py index a6ebd9498..3ad171429 100644 --- a/mesh/__init__.py +++ b/mesh/__init__.py @@ -3,4 +3,4 @@ necessary to work with finite-volume data. """ -__all__ = ['patch', 'integration', 'reconstruction'] +__all__ = ['integration', 'mapped', 'patch', 'reconstruction'] diff --git a/mesh/array_indexer.py b/mesh/array_indexer.py index 200c04c94..5dd7a6661 100644 --- a/mesh/array_indexer.py +++ b/mesh/array_indexer.py @@ -46,14 +46,14 @@ def __array_finalize__(self, obj): def __array_wrap__(self, out_arr, context=None): return np.ndarray.__array_wrap__(self, out_arr, context) - def v(self, buf=0, n=0, s=1): + def v(self, buf=0, n=None, s=1): """return a view of the valid data region for component n, with stride s, and a buffer of ghost cells given by buf """ return self.ip_jp(0, 0, buf=buf, n=n, s=s) - def ip(self, shift, buf=0, n=0, s=1): + def ip(self, shift, buf=0, n=None, s=1): """return a view of the data shifted by shift in the x direction. By default the view is the same size as the valid region, but the buf can specify how many ghost cells on each side to include. @@ -62,7 +62,7 @@ def ip(self, shift, buf=0, n=0, s=1): """ return self.ip_jp(shift, 0, buf=buf, n=n, s=s) - def jp(self, shift, buf=0, n=0, s=1): + def jp(self, shift, buf=0, n=None, s=1): """return a view of the data shifted by shift in the y direction. By default the view is the same size as the valid region, but the buf can specify how many ghost cells on each side to include. @@ -71,7 +71,7 @@ def jp(self, shift, buf=0, n=0, s=1): """ return self.ip_jp(0, shift, buf=buf, n=n, s=s) - def ip_jp(self, ishift, jshift, buf=0, n=0, s=1): + def ip_jp(self, ishift, jshift, buf=0, n=None, s=1): """return a view of the data shifted by ishift in the x direction and jshift in the y direction. By default the view is the same size as the valid region, but the buf can specify how many @@ -86,8 +86,12 @@ def ip_jp(self, ishift, jshift, buf=0, n=0, s=1): return np.asarray(self[self.g.ilo-bxlo+ishift:self.g.ihi+1+bxhi+ishift:s, self.g.jlo-bylo+jshift:self.g.jhi+1+byhi+jshift:s]) else: - return np.asarray(self[self.g.ilo-bxlo+ishift:self.g.ihi+1+bxhi+ishift:s, - self.g.jlo-bylo+jshift:self.g.jhi+1+byhi+jshift:s, n]) + if n is not None: + return np.asarray(self[self.g.ilo-bxlo+ishift:self.g.ihi+1+bxhi+ishift:s, + self.g.jlo-bylo+jshift:self.g.jhi+1+byhi+jshift:s, n]) + else: + return np.asarray(self[self.g.ilo-bxlo+ishift:self.g.ihi+1+bxhi+ishift:s, + self.g.jlo-bylo+jshift:self.g.jhi+1+byhi+jshift:s]) def lap(self, n=0, buf=0): """return the 5-point Laplacian""" diff --git a/mesh/integration.py b/mesh/integration.py index 7448c1c93..682560eef 100644 --- a/mesh/integration.py +++ b/mesh/integration.py @@ -101,13 +101,13 @@ def store_increment(self, istage, k_stage): weighting""" self.k[istage] = k_stage - def get_stage_start(self, istage): + def get_stage_start(self, istage, clone_function=patch.cell_center_data_clone): """get the starting conditions (a CellCenterData2d object) for stage istage""" if istage == 0: ytmp = self.start else: - ytmp = patch.cell_center_data_clone(self.start) + ytmp = clone_function(self.start) for n in range(ytmp.nvar): var = ytmp.get_var_by_index(n) for s in range(istage): diff --git a/mesh/mapped.py b/mesh/mapped.py new file mode 100644 index 000000000..a999938c9 --- /dev/null +++ b/mesh/mapped.py @@ -0,0 +1,434 @@ +""" +The patch module defines the classes necessary to describe finite-volume +data and the grid that it lives on. + +Typical usage: + +* create the grid:: + + grid = Grid2d(nx, ny) + +* create the data that lives on that grid:: + + data = CellCenterData2d(grid) + + bc = BC(xlb="reflect", xrb="reflect", + ylb="outflow", yrb="outflow") + data.register_var("density", bc) + ... + + data.create() + +* initialize some data:: + + dens = data.get_var("density") + dens[:, :] = ... + + +* fill the ghost cells:: + + data.fill_BC("density") + +""" +from __future__ import print_function + +import numpy as np +import sympy +from sympy.abc import x, y, z +from random import random +from numpy.testing import assert_array_almost_equal +from numba import njit + +from util import msg +import mesh.boundary as bnd +import mesh.array_indexer as ai +from mesh.patch import Grid2d, CellCenterData2d + + +class MappedGrid2d(Grid2d): + """ + the mapped 2-d grid class. The grid object will contain the coordinate + information (at various centerings). + + A basic (1-d) representation of the layout is:: + + | | | X | | | | X | | | + +--*--+- // -+--*--X--*--+--*--+- // -+--*--+--*--X--*--+- // -+--*--+ + 0 ng-1 ng ng+1 ... ng+nx-1 ng+nx 2ng+nx-1 + + ilo ihi + + |<- ng guardcells->|<---- nx interior zones ----->|<- ng guardcells->| + + The '*' marks the data locations. + """ + + def __init__(self, map_func, nx, ny, ng=1, + xmin=0.0, xmax=1.0, ymin=0.0, ymax=1.0): + """ + Create a MappedGrid2d object. + + The only data that we require is the number of points that + make up the mesh in each direction. Optionally we take the + extrema of the domain (default is [0,1]x[0,1]) and number of + ghost cells (default is 1). + + Note that the Grid2d object only defines the discretization, + it does not know about the boundary conditions, as these can + vary depending on the variable. + + Parameters + ---------- + map_func : sympy.Matrix + nx : int + Number of zones in the x-direction + ny : int + Number of zones in the y-direction + ng : int, optional + Number of ghost cells + xmin : float, optional + Mapped coordinate at the lower x boundary + xmax : float, optional + Mapped coordinate at the upper x boundary + ymin : float, optional + Mapped coordinate at the lower y boundary + ymax : float, optional + Mapped coordinate at the upper y boundary + """ + + super().__init__(nx, ny, ng, xmin, xmax, ymin, ymax) + + # we need to add a z-direction so that we can calculate the cross product + # of the basis vectors + self.map = map_func(self).col_join(sympy.Matrix([z])) + + print("Calculating scaling factors...") + self.kappa, self.gamma_fcx, self.gamma_fcy = self.calculate_metric_elements() + + print("Calculating rotation matrices...") + self.R_fcx, self.R_fcy = self.calculate_rotation_matrices() + + @staticmethod + def norm(z): + return sympy.sqrt(z.dot(z)) + + def sym_rotation_matrix(self): + """ + Use sympy to calculate the rotation matrices + """ + + Rx = sympy.zeros(2) + Ry = sympy.zeros(2) + + Rx[0, 0] = sympy.simplify(self.map[1].diff(y)) + Rx[0, 1] = sympy.simplify(self.map[0].diff(y)) + Rx[1, 0] = -sympy.simplify(self.map[0].diff(y)) + Rx[1, 1] = sympy.simplify(self.map[1].diff(y)) + + Ry[0, 0] = sympy.simplify(self.map[0].diff(x)) + Ry[0, 1] = sympy.simplify(self.map[1].diff(x)) + Ry[1, 0] = -sympy.simplify(self.map[1].diff(x)) + Ry[1, 1] = sympy.simplify(self.map[0].diff(x)) + + # normalize + Rx[0, :] /= self.norm(Rx[0, :]) + Rx[1, :] /= self.norm(Rx[1, :]) + Ry[0, :] /= self.norm(Ry[0, :]) + Ry[1, :] /= self.norm(Ry[1, :]) + + Rx = sympy.simplify(Rx) + Ry = sympy.simplify(Ry) + + # check rotation matrices - do this by substituting in random (non-zero) + # numbers as sympy is not great at cancelling things + assert_array_almost_equal((Rx @ Rx.T).subs( + {x: random() + 0.01, y: random() + 0.01}), np.eye(2)) + assert_array_almost_equal((Ry @ Ry.T).subs( + {x: random() + 0.01, y: random() + 0.01}), np.eye(2)) + + return sympy.simplify(Rx), sympy.simplify(Ry) + + def calculate_metric_elements(self): + """ + Given the functions for the area and line elements, calculate them on + the grid. + """ + + kappa = self.scratch_array() + hx = self.scratch_array() + hy = self.scratch_array() + + # calculate physical coordinates of corners + c1 = self.physical_coords(self.xl, self.yl) + c2 = self.physical_coords(self.xl, self.yr) + c3 = self.physical_coords(self.xr, self.yr) + c4 = self.physical_coords(self.xr, self.yl) + + def mapped_distance(m1, m2): + return np.sqrt((m1[0] - m2[0])**2 + (m1[1] - m2[1])**2) + + def mapped_area(i,j): + # find vectors of diagonals (and pad out z-direction with a zero) + p = np.append(c3[:,i,j] - c1[:,i,j], 0) + q = np.append(c4[:,i,j] - c2[:,i,j], 0) + + # area is half the cross product + return 0.5 * np.abs(np.cross(p, q)[-1]) + + for i in range(self.qx): + for j in range(self.qy): + hx[i, j] = mapped_distance(c1[:,i,j], c2[:,i,j]) + hy[i, j] = mapped_distance(c1[:,i,j], c4[:,i,j]) + kappa[i, j] = mapped_area(i, j) + + return kappa / (self.dx * self.dy), hx / self.dy, hy / self.dx + + def calculate_rotation_matrices(self): + """ + Calculate the rotation matrices on the cell interfaces. + It will return this as functions of nvar, ixmom and iymom - the grid + itself knows nothing of the variables, so these must be specified + by the MappedCellCenterData2d object. + """ + + # if isinstance(self.map, sympy.Matrix): + sym_Rx, sym_Ry = self.sym_rotation_matrix() + print('Rx = ', sym_Rx) + print('Ry = ', sym_Ry) + + Rx = sympy.lambdify((x, y), sym_Rx, modules="sympy") + Ry = sympy.lambdify((x, y), sym_Ry, modules="sympy") + + def R_fcx(nvar, ixmom, iymom): + R_fc = self.scratch_array(nvar=(nvar, nvar)) + + R_mat = np.eye(nvar) + + for i in range(self.qx): + for j in range(self.qy): + R_fc[i, j, :, :] = R_mat + + R_fc[i, j, ixmom:iymom + 1, ixmom:iymom + + 1] = np.array(Rx(self.xl[i], self.y[j])) + + return R_fc + + def R_fcy(nvar, ixmom, iymom): + R_fc = self.scratch_array(nvar=(nvar, nvar)) + + R_mat = np.eye(nvar) + + for i in range(self.qx): + for j in range(self.qy): + R_fc[i, j, :, :] = R_mat + + R_fc[i, j, ixmom:iymom + 1, ixmom:iymom + + 1] = np.array(Ry(self.x[i], self.yl[j])) + return R_fc + + return R_fcx, R_fcy + + def scratch_array(self, nvar=1): + """ + return a standard numpy array dimensioned to have the size + and number of ghostcells as the parent grid. + + Here I've generalized the version in Grid2d so that we can define + tensors (not just scalars and vectors) e.g. the rotation matrices + """ + + def flatten(t): + if not isinstance(t, tuple): + return (t, ) + elif len(t) == 0: + return () + else: + return flatten(t[0]) + flatten(t[1:]) + + if nvar == 1: + _tmp = np.zeros((self.qx, self.qy), dtype=np.float64) + else: + _tmp = np.zeros((self.qx, self.qy) + + flatten(nvar), dtype=np.float64) + return ai.ArrayIndexer(d=_tmp, grid=self) + + def physical_coords(self, xx=None, yy=None): + + if xx is None: + xs = self.x2d + elif not np.isscalar(xx) and len(xx.shape) == 1: + xs = np.repeat(xx, len(yy)) + xs.shape = (len(xx), len(yy)) + else: + xs = xx + + if yy is None: + ys = self.y2d + elif not np.isscalar(yy) and len(yy.shape) == 1: + ys = np.repeat(yy, len(xx)) + ys.shape = (len(yy), len(xx)) + ys = np.transpose(ys) + else: + ys = yy + + xs_t = sympy.lambdify((x, y), self.map[0]) + ys_t = sympy.lambdify((x, y), self.map[1]) + + return np.array([xs_t(xs, ys), ys_t(xs, ys)]) + + +class MappedCellCenterData2d(CellCenterData2d): + + def __init__(self, grid, dtype=np.float64): + + super().__init__(grid, dtype=dtype) + + self.R_fcx = [] + self.R_fcy = [] + + def make_rotation_matrices(self, ivars): + """ + The grid knows nothing of the variables, so we're going to define + the actual rotation matrices here by passing in the variable data + to the rotation matrix function. + """ + + self.R_fcx = self.grid.R_fcx(ivars.nvar, ivars.ixmom, ivars.iymom) + self.R_fcy = self.grid.R_fcy(ivars.nvar, ivars.ixmom, ivars.iymom) + + def fill_BC(self, name): + n = self.names.index(name) + self.fill_mapped_ghost(name, n=n, bc=self.BCs[name]) + + # that will handle the standard type of BCs, but if we asked + # for a custom BC, we handle it here + if self.BCs[name].xlb in bnd.ext_bcs.keys(): + bnd.ext_bcs[self.BCs[name].xlb](self.BCs[name].xlb, "xlb", name, self) + if self.BCs[name].xrb in bnd.ext_bcs.keys(): + bnd.ext_bcs[self.BCs[name].xrb](self.BCs[name].xrb, "xrb", name, self) + if self.BCs[name].ylb in bnd.ext_bcs.keys(): + bnd.ext_bcs[self.BCs[name].ylb](self.BCs[name].ylb, "ylb", name, self) + if self.BCs[name].yrb in bnd.ext_bcs.keys(): + bnd.ext_bcs[self.BCs[name].yrb](self.BCs[name].yrb, "yrb", name, self) + + def fill_mapped_ghost(self, name, n=0, bc=None): + """ + This replaces the fill_ghost routine in array_indexer to add some + scaling factors for mapped grids. + + We'll call fill_ghost then go back and fix the reflect-odd boundaries. + + It would be great if we could pass the variables object in here as we + actually only need to do this rotation if n == ixmom or iymom. + """ + myd = self.data + + myd.fill_ghost(n=n, bc=bc) + + if name == 'x-momentum' or name == 'y-momentum': + + # -x boundary + if bc.xlb in ["reflect-odd", "dirichlet"]: + if bc.xl_value is None: + for i in range(myd.g.ilo): + for j in range(myd.g.qy): + q_rot = self.R_fcy[myd.g.ilo, j] @ myd[2*myd.g.ng-i-1, j, :] + q_rot[n] = -q_rot[n] + myd[i, j, n] = (self.R_fcy[myd.g.ilo, j].T @ q_rot)[n] + else: + for j in range(myd.g.qy): + q_rot = self.R_fcy[myd.g.ilo, j] @ myd[myd.g.ilo, j, :] + q_rot[n] = 2*bc.xl_value - q_rot[n] + myd[myd.g.ilo-1, j, n] = (self.R_fcy[myd.g.ilo, j].T @ q_rot)[n] + + # +x boundary + if bc.xrb in ["reflect-odd", "dirichlet"]: + if bc.xr_value is None: + for i in range(myd.g.ng): + for j in range(myd.g.qy): + i_bnd = myd.g.ihi+1+i + i_src = myd.g.ihi-i + + q_rot = self.R_fcy[myd.g.ihi+1, j] @ myd[i_src, j, :] + q_rot[n] = -q_rot[n] + + myd[i_bnd, j, n] = (self.R_fcy[myd.g.ihi+1, j].T @ q_rot)[n] + else: + for j in range(myd.g.qy): + q_rot = self.R_fcy[myd.g.ihi+1, j] @ myd[myd.g.ihi, j, :] + q_rot[n] = 2*bc.xr_value - q_rot[n] + + myd[myd.g.ihi+1, j, n] = (self.R_fcy[myd.g.ihi+1, j].T @ q_rot)[n] + + # -y boundary + if bc.ylb in ["reflect-odd", "dirichlet"]: + if bc.yl_value is None: + for i in range(myd.g.qx): + for j in range(myd.g.jlo): + q_rot = self.R_fcx[i, myd.g.jlo] @ myd[i, 2*myd.g.ng-j-1, :] + q_rot[n] = -q_rot[n] + myd[i, j, n] = (self.R_fcx[i, myd.g.jlo].T @ q_rot)[n] + else: + for i in range(myd.g.qx): + q_rot = self.R_fcx[i, myd.g.jlo] @ myd[i, myd.g.jlo, :] + q_rot[n] = 2*bc.yl_value - q_rot[n] + myd[i, myd.g.jlo-1, n] = (self.R_fcx[i, myd.g.jlo].T @ q_rot)[n] + + # +y boundary + if bc.yrb in ["reflect-odd", "dirichlet"]: + if bc.yr_value is None: + for i in range(myd.g.qx): + for j in range(myd.g.ng): + j_bnd = myd.g.jhi+1+j + j_src = myd.g.jhi-j + + q_rot = self.R_fcx[i, myd.g.jhi+1] @ myd[i, j_src, :] + q_rot[n] = -q_rot[n] + + myd[i, j_bnd, n] = (self.R_fcx[i, myd.g.jhi+1].T @ q_rot)[n] + else: + for i in range(myd.g.qx): + + q_rot = self.R_fcx[i, myd.g.jhi+1] @ myd[i, myd.g.jhi, :] + q_rot[n] = 2*bc.yr_value - q_rot[n] + + myd[:, myd.g.jhi+1, n] = (self.R_fcx[i, myd.g.jhi+1].T @ q_rot)[n] + + +def mapped_cell_center_data_clone(old): + """ + Create a new CellCenterData2d object that is a copy of an existing + one + + Parameters + ---------- + old : CellCenterData2d object + The CellCenterData2d object we wish to copy + + Note + ---- + It may be that this whole thing can be replaced with a copy.deepcopy() + + """ + + if not isinstance(old, MappedCellCenterData2d): + msg.fail("Can't clone object") + + # we may be a type derived from CellCenterData2d, so use the same + # type + myt = type(old) + new = myt(old.grid, dtype=old.dtype) + + for n in range(old.nvar): + new.register_var(old.names[n], old.BCs[old.names[n]]) + + new.create() + + new.aux = old.aux.copy() + new.data = old.data.copy() + new.derives = old.derives.copy() + + new.R_fcx = old.R_fcx.copy() + new.R_fcy = old.R_fcy.copy() + + return new diff --git a/mesh/tests/test_mapped.py b/mesh/tests/test_mapped.py new file mode 100644 index 000000000..a6606a27e --- /dev/null +++ b/mesh/tests/test_mapped.py @@ -0,0 +1,174 @@ +# unit tests for mapped +import numpy as np +from numpy.testing import assert_array_equal, assert_array_almost_equal +import sympy +from sympy.abc import x, y + +import mesh.mapped as mapped + + +class TestMappedGrid2d(object): + + @classmethod + def setup_class(cls): + """ this is run once for each class before any tests """ + pass + + @classmethod + def teardown_class(cls): + """ this is run once for each class after all tests """ + pass + + def setup_method(self): + """ this is run before each test """ + pass + + def teardown_method(self): + """ this is run after each test """ + pass + + def test_rectilinear(self): + """ + Test mapped grid class for a rectilinear grid + """ + + def map(myg): + return sympy.Matrix([2 * x, 0.1 * y]) + + grid = mapped.MappedGrid2d(map, 4, 8, ng=2, ymax=2.0) + + assert grid.dx == 0.25 + assert grid.dy == 0.25 + + assert_array_almost_equal( + grid.kappa, np.ones((8, 12)) * 0.2, decimal=12) + assert_array_almost_equal( + grid.gamma_fcx, np.ones((8, 12)) * 0.1, decimal=12) + assert_array_almost_equal( + grid.gamma_fcy, np.ones((8, 12)) * 2, decimal=12) + + assert_array_almost_equal(grid.R_fcx( + 2, 0, 1)[0, 0, :, :], np.eye(2), decimal=12) + assert_array_almost_equal(grid.R_fcy( + 2, 0, 1)[-1, -1, :, :], np.eye(2), decimal=12) + + def test_noncartesian(self): + """ + Test mapped grid class for a non cartesian grid. + """ + + def map(myg): + return sympy.Matrix([0.5 * x + y, x - 0.4 * y]) + + grid = mapped.MappedGrid2d(map, 4, 4, ng=2) + + # numpy and scipy use slightly different sqrt routines so lose a little + # bit of precision - use almost equal here rather than equal + assert_array_almost_equal(grid.kappa, np.ones((8, 8)) * 1.2) + assert_array_almost_equal(grid.gamma_fcx, np.ones( + (8, 8)) * np.sqrt(0.4**2 + 1), decimal=12) + assert_array_almost_equal(grid.gamma_fcy, np.ones( + (8, 8)) * np.sqrt(0.5**2 + 1), decimal=12) + + Rx = np.array([[-0.4 / np.sqrt(0.4**2 + 1), 1 / np.sqrt(0.4**2 + 1)], + [-1 / np.sqrt(0.4**2 + 1), -0.4 / np.sqrt(0.4**2 + 1)]]) + + Ry = np.array([[0.5 / np.sqrt(0.5**2 + 1), 1 / np.sqrt(0.5**2 + 1)], + [-1 / np.sqrt(0.5**2 + 1), 0.5 / np.sqrt(0.5**2 + 1)]]) + + assert_array_almost_equal(grid.R_fcx( + 2, 0, 1)[0, 0, :, :], Rx, decimal=12) + assert_array_almost_equal(grid.R_fcy( + 2, 0, 1)[2, 4, :, :], Ry, decimal=12) + + def test_polar(self): + """ + Test mapped grid class for polar coordinates. + """ + + def map(myg): + return sympy.Matrix([x * sympy.cos(y), x * sympy.sin(y)]) + + grid = mapped.MappedGrid2d(map, 4, 4, ng=2, xmin=0.1) + + xs = grid.x2d * np.cos(grid.y2d) + ys = grid.x2d * np.sin(grid.y2d) + + x_map = sympy.lambdify((x, y), map(grid)[0])(grid.x2d, grid.y2d) + y_map = sympy.lambdify((x, y), map(grid)[1])(grid.x2d, grid.y2d) + + assert_array_almost_equal(xs, x_map, decimal=12) + assert_array_almost_equal(ys, y_map, decimal=12) + + kappa = np.array([[0.23503376, 0.23503376, 0.23503376, 0.23503376, 0.23503376, 0.23503376, 0.23503376, 0.23503376], + [0.0123702, 0.0123702, 0.0123702, 0.0123702, 0.0123702, 0.0123702, + 0.0123702, 0.0123702], + [0.21029337, 0.21029337, 0.21029337, 0.21029337, 0.21029337, 0.21029337, + 0.21029337, 0.21029337], + [0.43295693, 0.43295693, 0.43295693, 0.43295693, 0.43295693, 0.43295693, + 0.43295693, 0.43295693], + [0.65562049, 0.65562049, 0.65562049, 0.65562049, 0.65562049, 0.65562049, + 0.65562049, 0.65562049], + [0.87828406, 0.87828406, 0.87828406, 0.87828406, 0.87828406, 0.87828406, + 0.87828406, 0.87828406], + [1.10094762, 1.10094762, 1.10094762, 1.10094762, 1.10094762, 1.10094762, + 1.10094762, 1.10094762], + [1.32361118, 1.32361118, 1.32361118, 1.32361118, 1.32361118, 1.32361118, + 1.32361118, 1.32361118]]) + assert_array_almost_equal(grid.kappa, kappa) + + gamma_fcx = np.array([[0.34908925, 0.34908925, 0.34908925, 0.34908925, 0.34908925, 0.34908925, 0.34908925, 0.34908925], + [0.12467473, 0.12467473, 0.12467473, 0.12467473, + 0.12467473, 0.12467473, 0.12467473, 0.12467473], + [0.09973979, 0.09973979, 0.09973979, 0.09973979, 0.09973979, 0.09973979, + 0.09973979, 0.09973979], + [0.32415431, 0.32415431, 0.32415431, 0.32415431, 0.32415431, 0.32415431, + 0.32415431, 0.32415431], + [0.54856883, 0.54856883, 0.54856883, 0.54856883, 0.54856883, 0.54856883, + 0.54856883, 0.54856883], + [0.77298335, 0.77298335, 0.77298335, 0.77298335, 0.77298335, 0.77298335, + 0.77298335, 0.77298335], + [0.99739787, 0.99739787, 0.99739787, 0.99739787, 0.99739787, 0.99739787, + 0.99739787, 0.99739787], + [1.22181239, 1.22181239, 1.22181239, 1.22181239, 1.22181239, 1.22181239, + 1.22181239, 1.22181239]] + ) + assert_array_almost_equal(grid.gamma_fcx, gamma_fcx) + + gamma_fcy = np.ones_like(xs) + assert_array_almost_equal(grid.gamma_fcy, gamma_fcy, decimal=12) + + Rx = lambda xx, yy: np.array([[np.cos(yy), -np.sin(yy)], + [np.sin(yy), np.cos(yy)]]) + + Ry = lambda xx, yy: np.array([[np.cos(yy), np.sin(yy)], + [-np.sin(yy), np.cos(yy)]]) + + R_fcx = Rx(grid.x2d[4, 5] - 0.5 * grid.dx, grid.y2d[4, 5]) + R_fcy = Ry(grid.x2d[3, 2], grid.y2d[3, 2] - 0.5 * grid.dy) + + assert_array_almost_equal(grid.R_fcx( + 2, 0, 1)[4, 5, :, :], R_fcx, decimal=12) + assert_array_almost_equal(grid.R_fcy( + 2, 0, 1)[3, 2, :, :], R_fcy, decimal=12) + + +class TestMappedCellCenterData2d(object): + + @classmethod + def setup_class(cls): + """ this is run once for each class before any tests """ + pass + + @classmethod + def teardown_class(cls): + """ this is run once for each class after all tests """ + pass + + def setup_method(self): + """ this is run before each test """ + self.g = mapped.MappedGrid2d(4, 6, ng=2, ymax=1.5) + + def teardown_method(self): + """ this is run after each test """ + self.g = None diff --git a/pyro.py b/pyro.py index 67de7b688..a2dcc1657 100755 --- a/pyro.py +++ b/pyro.py @@ -19,6 +19,7 @@ "compressible", "compressible_rk", "compressible_fv4", + "compressible_mapped", "compressible_sdc", "compressible_react", "diffusion", diff --git a/simulation_null.py b/simulation_null.py index 1f98ca39d..5bc80a941 100644 --- a/simulation_null.py +++ b/simulation_null.py @@ -82,6 +82,20 @@ def bc_setup(rp): return bc, bc_xodd, bc_yodd +class NullVariables(object): + """ + a container class for easy access to the different variables by an integer key + """ + + def __init__(self, myd): + self.nvar = len(myd.names) + + self.initialize(myd) + + def initialize(self, myd): + pass + + class NullSimulation(object): def __init__(self, solver_name, problem_name, rp, timers=None, data_class=patch.CellCenterData2d): @@ -153,7 +167,8 @@ def do_output(self): n_out = self.rp.get_param("io.n_out") do_io = self.rp.get_param("io.do_io") - is_time = self.cc_data.t >= (self.n_num_out + 1)*dt_out or self.n % n_out == 0 + is_time = self.cc_data.t >= ( + self.n_num_out + 1) * dt_out or self.n % n_out == 0 if is_time and do_io == 1: self.n_num_out += 1 return True @@ -185,9 +200,9 @@ def compute_timestep(self): else: self.method_compute_timestep() if self.n == 0: - self.dt = init_tstep_factor*self.dt + self.dt = init_tstep_factor * self.dt else: - self.dt = min(max_dt_change*self.dt_old, self.dt) + self.dt = min(max_dt_change * self.dt_old, self.dt) self.dt_old = self.dt if self.cc_data.t + self.dt > self.tmax: @@ -215,7 +230,8 @@ def finalize(self): finalize() method. """ # there should be a cleaner way of doing this - problem = importlib.import_module("{}.problems.{}".format(self.solver_name, self.problem_name)) + problem = importlib.import_module( + "{}.problems.{}".format(self.solver_name, self.problem_name)) problem.finalize()