rotation matrix is less broken

compressible/problems/acoustic_pulse.py

@@ -12,11 +12,11 @@ def init_data(myd, rp):
 
     msg.bold("initializing the acoustic pulse problem...")
 
-    # # make sure that we are passed a valid patch object
-    # if not isinstance(myd, fv.FV2d):
-    #     print("ERROR: patch invalid in acoustic_pulse.py")
-    #     print(myd.__class__)
-    #     sys.exit()
+    # make sure that we are passed a valid patch object
+    if not isinstance(myd, fv.FV2d):
+        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")

compressible_mapped/problems/acoustic_pulse.py

@@ -5,7 +5,7 @@
 import numpy as np
 from util import msg
 import sympy
-from sympy.abc import x, y
+from sympy.abc import x, y, p, q
 
 
 def init_data(myd, rp):
@@ -45,11 +45,14 @@ def init_data(myd, rp):
 
     myg = myd.grid
 
-    xctr = 0.5 * (xmin + xmax) * myg.gamma_fcy
-    yctr = 0.5 * (ymin + ymax) * myg.gamma_fcx
+    xctr = 0.5 * (xmin + xmax)
+    yctr = 0.5 * (ymin + ymax)
 
-    dist = np.sqrt((myd.grid.x2d * myg.gamma_fcy - xctr)**2 +
-                   (myd.grid.y2d * myg.gamma_fcx - yctr)**2)
+    X, Y = myg.physical_coords()
+
+    xctr_t, yctr_t = myg.physical_coords(xctr, yctr)
+
+    dist = np.sqrt((X - xctr_t)**2 + (Y - yctr_t)**2)
 
     dens[:, :] = rho0
     idx = dist <= 0.5
@@ -60,40 +63,13 @@ def init_data(myd, rp):
     ener[:, :] = p / (gamma - 1)
 
 
-def area(myg):
-    return 2 * myg.dx * myg.dy + myg.scratch_array()
-
-
-def h(idir, myg):
-    if idir == 1:
-        return myg.dy + myg.scratch_array()
-    else:
-        return 2 * myg.dx + myg.scratch_array()
-
-
-def R(iface, myg, nvar, ixmom, iymom):
-    R_fc = myg.scratch_array(nvar=(nvar, nvar))
-
-    R_mat = np.eye(nvar)
-
-    for i in range(myg.qx):
-        for j in range(myg.qy):
-            R_fc[i, j, :, :] = R_mat
-
-    return R_fc
-
-
-def map(myg):
-
-    xs_t = myg.x2d * 2
-    ys_t = myg.y2d
-
-    return xs_t, ys_t
+def sym_map(myg):
 
+    return sympy.Matrix([x+3 + y, -2*y])
 
-def sym_map(myg):
+def inverse_map(myg):
 
-    return sympy.Matrix([2 * x, y-1])
+    return sympy.Matrix([])
 
 
 def finalize():

compressible_mapped/simulation.py

@@ -217,16 +217,13 @@ def dovis(self):
 
         _, 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]
 
-            X = myg.x2d[myg.ng:-myg.ng, myg.ng:-myg.ng] * myg.gamma_fcy.v()
-            Y = myg.y2d[myg.ng:-myg.ng, myg.ng:-myg.ng] * myg.gamma_fcx.v()
-
-            # X, Y = myg.map
-            # X = X[myg.ng:-myg.ng, myg.ng:-myg.ng]
-            # Y = Y[myg.ng:-myg.ng, myg.ng:-myg.ng]
-
             img = ax.pcolormesh(X, Y, v.v(), cmap=self.cm)
 
             ax.set_xlabel("x")

mesh/mapped.py

@@ -137,16 +137,27 @@ def sym_line_elements(self):
 
     def sym_rotation_matrix(self):
         """
-        Use sympy to calculate the rotation matrix
+        Use sympy to calculate the rotation matrices
         """
 
         J = sympy.Matrix(self.map[:-1]).jacobian((x, y))
-        if J[0, :].norm() != 0:
-            J[0, :] /= J[0, :].norm()
-        if J[1, :].norm() != 0:
-            J[1, :] /= J[1, :].norm()
 
-        return J
+        # normalize
+        J[0,:] /= J[0,:].norm()
+        J[1,:] /= J[1,:].norm()
+
+        Rx = sympy.zeros(2)
+        Ry = sympy.zeros(2)
+
+        Rx[1,:] = J[1,:]
+        Rx[0,0] = J[1,1]
+        Rx[0,1] = -J[1,0]
+
+        Ry[1,:] = J[0,::-1]
+        Ry[0,0] = J[0,0]
+        Ry[0,1] = -J[0,1]
+
+        return Rx, Ry
 
     def calculate_metric_elements(self):
         """
@@ -180,6 +191,10 @@ def calculate_metric_elements(self):
         #     hx[:, :] = h(1, self) / self.dy
         #     hy[:, :] = h(2, self) / self.dx
 
+        print('dA = ', sym_dA)
+        print('hx = ', sym_hx)
+        print('hy = ', sym_hy)
+
         return kappa, hx, hy
 
     def calculate_rotation_matrices(self):
@@ -191,45 +206,47 @@ def calculate_rotation_matrices(self):
         """
 
         # if isinstance(self.map, sympy.Matrix):
-        sym_R = self.sym_rotation_matrix()
-        R = sympy.lambdify((x, y), sym_R)
+        sym_Rx, sym_Ry = self.sym_rotation_matrix()
+        print('Rx = ', sym_Rx)
+        print('Ry = ', sym_Ry)
+
+        # R = sympy.lambdify((x, y), sym_R, modules="sympy")
+
+        # print(sympy.limit(sym_R, x, 0))
 
         def R_fcx(nvar, ixmom, iymom):
             R_fc = self.scratch_array(nvar=(nvar, nvar))
 
             R_mat = np.eye(nvar)
 
+            xs = self.x2d - 0.5 * self.dx
+            ys = self.y2d
+
             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] = R(self.x2d[i, j] - 0.5 * self.dx, self.y2d[i, j])
+                         1] = np.array(sym_Rx.subs({x: xs[i, j], y: ys[i, j]}))
 
             return R_fc
 
-        # print(R_fcx(4, 2, 3))
-
         def R_fcy(nvar, ixmom, iymom):
             R_fc = self.scratch_array(nvar=(nvar, nvar))
 
             R_mat = np.eye(nvar)
 
+            xs = self.x2d
+            ys = self.y2d - 0.5 * self.dy
+
             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] = R(self.x2d[i, j], self.y2d[i, j] - 0.5 * self.dy)
+                         1] = np.array(sym_Ry.subs({x: xs[i, j], y: ys[i, j]}))
             return R_fc
 
-            # print(R_fcy(4, 2, 3))
-        # else:
-        #     R_fcx = partial(R, 1, self)
-        #     R_fcy = partial(R, 2, self)
-
-            # print(R_fcx(4, 2, 3))
-
         return R_fcx, R_fcy
 
     def scratch_array(self, nvar=1):
@@ -256,6 +273,18 @@ def flatten(t):
                             flatten(nvar), dtype=np.float64)
         return ai.ArrayIndexer(d=_tmp, grid=self)
 
+    def physical_coords(self, xs=None, ys=None):
+
+        if xs is None:
+            xs = self.x2d
+        if ys is None:
+            ys = self.y2d
+
+        xs_t = sympy.lambdify((x, y), self.map[0])
+        ys_t = sympy.lambdify((x, y), self.map[1])
+
+        return xs_t(xs, ys), ys_t(xs, ys)
+
 
 class MappedCellCenterData2d(CellCenterData2d):
 
@@ -276,6 +305,8 @@ def make_rotation_matrices(self, ivars):
         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)
 
+        # print('Rx contains nan?', np.isnan(self.R_fcx).any())
+
 
 def mapped_cell_center_data_clone(old):
     """