From ef42b3754ef83f5ba817b20a135a93e46af18d71 Mon Sep 17 00:00:00 2001 From: Cory Mikida Date: Tue, 31 Aug 2021 15:45:57 -0500 Subject: [PATCH 1/7] Towards SBP-DG interface (brought over from Gitlab) --- grudge/models/advection.py | 49 ++++ test/projection.py | 539 +++++++++++++++++++++++++++++++++++++ test/q_2.txt | 6 + test/q_4.txt | 96 +++++++ test/q_6.txt | 324 ++++++++++++++++++++++ test/sbp_operators.py | 163 +++++++++++ test/test_grudge.py | 368 +++++++++++++++++++++++++ test/test_sbp_dg.py | 393 +++++++++++++++++++++++++++ 8 files changed, 1938 insertions(+) create mode 100644 test/projection.py create mode 100644 test/q_2.txt create mode 100644 test/q_4.txt create mode 100644 test/q_6.txt create mode 100644 test/sbp_operators.py create mode 100644 test/test_sbp_dg.py diff --git a/grudge/models/advection.py b/grudge/models/advection.py index aacd2fcfd..10e7405d5 100644 --- a/grudge/models/advection.py +++ b/grudge/models/advection.py @@ -178,6 +178,55 @@ def flux(tpair): ) ) + +class WeakAdvectionSBPOperator(AdvectionOperatorBase): + def flux(self, u_tpair): + return self.weak_flux(u_tpair) + + from meshmode.mesh import BTAG_NONE + + def operator(self, t, u, state_from_sbp, + sbp_tag=BTAG_NONE, std_tag=BTAG_NONE): + + dcoll = self.dcoll + + def flux(tpair): + return op.project(dcoll, tpair.dd, "all_faces", self.flux(tpair)) + + if self.inflow_u is not None: + inflow_flux = flux(op.bv_trace_pair(dcoll, + std_tag, + interior=u, + exterior=self.inflow_u(t))) + else: + inflow_flux = 0 + + sbp_flux = flux(op.bv_trace_pair(dcoll, + sbp_tag, + interior=u, + exterior=state_from_sbp)) + return ( + op.inverse_mass( + dcoll, + np.dot(self.v, op.weak_local_grad(dcoll, u)) + - op.face_mass( + dcoll, + sum(flux(tpair) for tpair in op.interior_trace_pairs(dcoll, u)) + + inflow_flux + sbp_flux + + # FIXME: Add support for inflow/outflow tags + # + flux(op.bv_trace_pair(dcoll, + # self.inflow_tag, + # interior=u, + # exterior=bc_in)) + # + flux(op.bv_trace_pair(dcoll, + # self.outflow_tag, + # interior=u, + # exterior=bc_out)) + ) + ) + ) + # }}} diff --git a/test/projection.py b/test/projection.py new file mode 100644 index 000000000..c3e1b5a59 --- /dev/null +++ b/test/projection.py @@ -0,0 +1,539 @@ +import numpy as np +from sbp_operators import (sbp21, sbp42, sbp63) + + +def gaussian_quad(n): + + # Gaussian quadrature per the algorithm + # of Hesthaven and Warburton. + + if n == 0: + x = 0 + w = 2 + else: + h1 = 2 * np.linspace(0, n, n+1, endpoint=True) + je = 2/(h1[0:n] + 2)*np.linspace(1, n, n, endpoint=True) * \ + np.linspace(1, n, n, endpoint=True) / \ + np.sqrt((h1[0:n]+1)*(h1[0:n]+3)) + a = np.diag(je, 1) + ap = np.diag(je, -1) + a = a + ap + [x, v] = np.linalg.eig(a) + idx = x.argsort()[::1] + x = x[idx] + v = v[:, idx] + w = 2 * (v[0, :]) ** 2 + + return x, w + + +def legendre_vandermonde(x, n): + + v = np.zeros((x.shape[0], n+1)) + p_n_0 = x + p_n_1 = np.ones(x.shape[0]) + if n >= 0: + v[:, 0] = p_n_1 * np.sqrt(1.0/2.0) + + if n >= 1: + v[:, 1] = p_n_0 * np.sqrt(3.0/2.0) + + for i in range(2, n+1): + a = (2*i - 1)/i + c = (i - 1)*(i - 1)*(2*i) / (i * i * (2*i - 2)) + p_n_2 = p_n_1 + p_n_1 = p_n_0 + p_n_0 = a*x*p_n_1 - c*p_n_2 + v[:, i] = p_n_0 * np.sqrt((2*i+1)/2) + + return v + + +def glue_pieces(n, vx, xg): + + # Remove leading dimensions of 1. + if vx.shape[0] == 1: + vx = np.squeeze(vx, axis=0) + k = vx.shape[0] - 1 + + [r, w_dump] = gaussian_quad(n) + # Python uses zero-based indices. + va = np.arange(0, k, 1) + vb = np.arange(1, k+1, 1) + + # x = np.ones(n+1)*vx[va] + 0.5*(r+1)*(vx[vb] - vx[va]) + x = (np.ones(n+1).reshape(1, -1)).transpose()*vx[va] + \ + 0.5*((r.reshape(1, -1)).transpose()+1)*(vx[vb] - vx[va]) + x = np.squeeze(x) + + m = np.diag(2/(2*np.arange(0, n+1, 1)+1)) + vr = legendre_vandermonde(r, n) + u = np.zeros((n+1, k*(n+1))) + + xmin = vx[0] + xmax = vx[-1] + + xbr = 2*np.divide((x-xmin), (xmax-xmin)) - 1 + vbr = (legendre_vandermonde((xbr.transpose()).flatten(), + n)).dot(np.sqrt(m)) + + invsqm_invvr = np.sqrt(np.diag( + np.divide(1, np.diag(m)))).dot(np.linalg.inv(vr)) + + for i_n in range(0, n+1): + pbr = np.reshape(vbr[:, i_n], x.shape) + tmp = invsqm_invvr.dot(pbr) + u[i_n, :] = (tmp.transpose()).flatten() + + xgbr = 2*np.divide((xg-xmin), (xmax-xmin)) - 1 + v = (legendre_vandermonde(xgbr, n)).dot(np.sqrt(m)) + + m = np.divide(m, 2) + + return u, v, m + + +def make_projection(n, order): + + p = order - 1 + + # Diagonal SBP operators + # Just orders 2 and 4 for now. + if order == 2: + # In K+W code, FD op size is n+1 + [p_sbp, q_sbp] = sbp21(n+1) + m = 2 + # pb = 1 + # from opt files + q = np.loadtxt('q_2.txt') + r = 1 + elif order == 4: + # In K+W code, FD op size is n+1 + [p_sbp, q_sbp] = sbp42(n+1) + m = 4 + # pb = 2 + # from opt files + q = np.loadtxt('q_4.txt') + r = 5 + elif order == 6: + # In K+W code, FD op size is n+1 + [p_sbp, q_sbp] = sbp63(n+1) + m = 6 + # pb = 2 + # from opt files + q = np.loadtxt('q_6.txt') + r = 8 + + h = p_sbp + + s = int(r - (m/2 - 1)) + + # Make glue and finite difference grids. + xf_b = np.linspace(0, s + m - 2, s + m - 1, endpoint=True) + xg_b = np.linspace(0, s + m - 2, s + m - 1, endpoint=True) + # Get the glue pieces. + [u_dump, v_dump, mr_b] = glue_pieces(p, xf_b, xg_b) + + i = np.kron(np.linspace(s+1, n-s+1, n-2*s+1, endpoint=True), + np.ones(m*(p+1))) + j = np.kron(np.ones(n+1-2*s), + np.linspace(1, m*(p+1), m*(p+1), endpoint=True)) + \ + (p+1)*(i-1-m/2) + # Interior solution. + qi = np.kron(np.ones(n+1-2*s), q[0:(m*(p+1))]) + pg2f = np.zeros((n+1, n*(p+1))) + for k in range(0, i.shape[0]): + # Zero-based indices. + pg2f[int(i[k])-1, int(j[k])-1] = qi[k] + + # Boundary solution. + qb = np.reshape(q[m*(p+1):], (r*(p+1), s,), order='F') + qb = np.transpose(qb) + + # Left block. + pg2f[0:s, 0:r*(p+1)] = qb + + qb = np.reshape((np.diag(2*np.mod(np.linspace(1, p+1, p+1, endpoint=True), + 2) - 1)).dot(np.flipud(np.reshape(np.rot90(qb, + 2).transpose(), (p+1, r*s,), order='F'))), + (r*(p+1), s,), order='F').transpose() + + # pg2f[np.arange(pg2f.shape[0]-s, pg2f.shape[0], 1), + # np.arange(pg2f.shape[1]+1-r*(p+1)-1, pg2f.shape[1], 1)] = qb + + # Right block. + for ind_i in range(pg2f.shape[0]-s, pg2f.shape[0]): + for ind_j in range(pg2f.shape[1]+1-r*(p+1)-1, pg2f.shape[1]): + pg2f[ind_i, ind_j] = qb[ind_i-(pg2f.shape[0]-s), + ind_j - (pg2f.shape[1]+1-r*(p+1)-1)] + + m = np.kron(np.eye(n), mr_b) + + # Pf2g comes directly from the compatibility equation. + pf2g = (np.kron(np.eye(n), + np.diag(1/np.diag(mr_b))).dot(np.transpose(pg2f))).dot(h) + + return pf2g, pg2f, m, h + + +def make_projection_g2g_hr(n, vxi_l, vxi_r): + + n_p = n + 1 + + nv_l = vxi_l.shape[0] + nv_r = vxi_r.shape[0] + + k_l = nv_l - 1 + k_r = nv_r - 1 + + tol = 100.0*np.finfo(float).eps + + # Check to ensure that first and last grid points align. + assert abs(vxi_l[-1] - vxi_r[-1]) <= \ + (tol*max(abs(vxi_l[-1]), abs(vxi_r[-1])) + np.finfo(float).eps) + assert abs(vxi_l[-1] - vxi_r[-1]) <= \ + (tol*max(abs(vxi_l[-1]), abs(vxi_r[-1])) + np.finfo(float).eps) + + vxi_g = np.sort(np.concatenate([vxi_l, vxi_r])) + vxi_g = np.unique(vxi_g) + + xi_g = np.vstack((vxi_g[0:-1], vxi_g[1:])) + k_g = xi_g.shape[1] + + np_l = k_l * n_p + np_r = k_r * n_p + np_g = k_g * n_p + + pg2l = np.zeros((np_l, np_g,)) + pl2g = np.zeros((np_g, np_l,)) + + for k in range(0, k_l): + xia_l = vxi_l[k] + xib_l = vxi_l[k+1] + + k_g = np.where(np.all( + np.vstack(([vxi_g < xib_l - tol], + [vxi_g >= xia_l - tol])), axis=0)) + + [pc2f, pf2c] = make_projection_g2g_h_gen(n, np.squeeze(xi_g[:, k_g])) + + idx_l = (k) * n_p + idx_g = (k_g[0][0]) * n_p + + # pg2l[np.arange(idx_l, idx_l+n_p, 1), + # np.arange(idx_g, idx_g+n_p*k_g[0].shape[0], 1)] = pf2c + # pl2g[np.arange(idx_g, idx_g+n_p*kg[0].shape[0], 1), + # np.arange(idx_l, idx_l+n_p, 1)] = pc2f + + for i in range(idx_l, idx_l+n_p): + for j in range(idx_g, idx_g+n_p*k_g[0].shape[0]): + pg2l[i, j] = pf2c[i-idx_l, j-idx_g] + pl2g[j, i] = pc2f[j-idx_g, i-idx_l] + + pg2r = np.zeros((np_r, np_g)) + pr2g = np.zeros((np_g, np_r)) + + for k in range(0, k_r): + xia_r = vxi_r[k] + xib_r = vxi_r[k+1] + + k_g = np.where(np.all( + np.vstack(([vxi_g < xib_r - tol], + [vxi_g >= xia_r - tol])), axis=0)) + + [pc2f, pf2c] = make_projection_g2g_h_gen(n, np.squeeze(xi_g[:, k_g])) + + idx_r = (k) * n_p + idx_g = (k_g[0][0]) * n_p + + # pg2r[np.arange(idx_r, idx_r+n_p-1, 1), + # np.arange(idx_g, idx_g+n_p*k_g[0].shape[0]-1, 1)] = pf2c + # pr2g[idx_g:idx_g+n_p*k_g[0].shape[0]-1, + # idx_r:idx_r+n_p-1] = pc2f + + for i in range(idx_r, idx_r+n_p): + for j in range(idx_g, idx_g+n_p*k_g[0].shape[0]): + pg2r[i, j] = pf2c[i-idx_r, j-idx_g] + pr2g[j, i] = pc2f[j-idx_g, i-idx_r] + + return vxi_g, pl2g, pg2l, pr2g, pg2r + + +def make_projection_g2g_h_gen(n, xi_c): + + [r, w_dump] = gaussian_quad(n) + + if len(xi_c.shape) == 1: + k = 1 + vx = xi_c + fa = xi_c[0] + fb = xi_c[-1] + else: + k = xi_c.shape[1] + vx = np.append(xi_c[0, :], xi_c[1, -1]) + fa = xi_c[0, 0] + fb = xi_c[1, -1] + + vx = 2/(fb - fa) * vx + 1 - 2*fb/(fb-fa) + + va = np.arange(0, k, 1) + vb = np.arange(1, k+1, 1) + + x = (np.ones(n+1).reshape(1, -1)).transpose()*vx[va] + \ + 0.5*((r.reshape(1, -1)).transpose()+1)*(vx[vb] - vx[va]) + + jc = (vx[vb] - vx[va])/2 + jf = 1 + + pc2f = np.zeros((x.size, r.shape[0])) + + m = (2.0/(2*np.linspace(0, n, n+1, endpoint=True) + 1)) + + vr = legendre_vandermonde(r, n) + + sq_invm_invvr = np.diag(np.sqrt(np.divide(1, m))).dot(np.linalg.inv(vr)) + + for k_i in range(0, k): + pc2f[np.arange(k_i*(n+1), (k_i+1)*(n+1), 1), :] = \ + sq_invm_invvr.dot(legendre_vandermonde(x[:, k_i], + n).dot(np.diag(np.sqrt(m)))) + + pf2c = (1.0 / jf) * (np.diag(1.0 / m).dot( + pc2f.transpose())).dot(np.kron(np.diag(jc), np.diag(m))) + + return pc2f, pf2c + + +def sbp_sbp_projection(na, qa, nb, qb): + + # Get the projection operators for each side. + [paf2g, pag2f, ma_dump, ha_dump] = make_projection(na, qa) + [pbf2g, pbg2f, mb_dump, hb_dump] = make_projection(nb, qb) + + # Modify the projection operators so they go to the + # same order polynomials. + + pg = max(qa, qb)-1 + + paf2g = np.kron(np.eye(na), np.eye(pg+1, qa)).dot(paf2g) + pag2f = pag2f.dot(np.kron(np.eye(na), np.eye(qa, pg+1))) + + pbf2g = np.kron(np.eye(nb), np.eye(pg+1, qb)).dot(pbf2g) + pbg2f = pbg2f.dot(np.kron(np.eye(nb), np.eye(qb, pg+1))) + + # Move to the glue space + xa = np.linspace(-1, 1, na + 1, endpoint=True) + xb = np.linspace(-1, 1, nb + 1, endpoint=True) + + # Glue-to-glue. + [vxi_g_dump, pa2g, pg2a, pb2g, pg2b] = make_projection_g2g_hr(pg, xa, xb) + + # Stack the operators. + pa2b = ((pbg2f.dot(pg2b)).dot(pa2g)).dot(paf2g) + pb2a = ((pag2f.dot(pg2a)).dot(pb2g)).dot(pbf2g) + + return pa2b, pb2a + + +def sbp_sbp_test(): + + na = 100 + nb = 230 + + xa = np.linspace(-1, 1, na+1, endpoint=True) + xb = np.linspace(-1, 1, nb+1, endpoint=True) + + eps = np.finfo(float).eps + + # FIXME: these test ranges are supposed to go to 10, but + # we don't have those diagonal SBP operators coded yet. + for qa in range(2, 6, 2): + for qb in range(2, 6, 2): + print('Creating projection for (qa, qb) ', qa, qb) + [pa2b, pb2a] = sbp_sbp_projection(na, qa, nb, qb) + + # np.savetxt('pa2b_test.txt', pa2b) + # np.savetxt('pb2a_test.txt', pb2a) + + print('Testing projection for (qa, qb) ', qa, qb) + # Check the full operator + for n in range(0, int(min(qa, qb)/2)-1): + np.testing.assert_array_less(np.abs(pa2b.dot((xa ** n)) + - xb ** n), 1000*eps) + np.testing.assert_array_less(np.abs(pb2a.dot((xb ** n)) + - xa ** n), 1000*eps) + print('Test passed for (qa, qb) ', qa, qb) + + # Check the interior operator + n_int = int(min(qa, qb))-1 + ta = pb2a.dot(xb ** n_int) - xa ** n_int + tb = pa2b.dot(xa ** n_int) - xb ** n_int + + mb = np.argwhere(np.abs(tb) > 1000*eps).size / 2 + ma = np.argwhere(np.abs(ta) > 1000*eps).size / 2 + + assert mb < nb / 2 + assert ma < na / 2 + + # Look at interior part only - locate boundary portions ad-hoc. + if mb > 0: + assert np.max(abs(tb[int(mb):int(-mb)])) < 1000*eps + if ma > 0: + assert np.max(abs(ta[int(ma):int(-ma)])) < 1000*eps + + return + + +def sbp_dg_projection(na, nb, qa, qb, gg, dg_nodes): + + # Get the projection operators for SBP-side. + [paf2g, pag2f, ma_dump, ha_dump] = make_projection(na, qa) + [pbf2g, pbg2f] = glue_to_dg(qb, gg, dg_nodes) + + # Modify the projection operator so that it goes to the + # same order polynomials as DG. + + pg = max(qa, qb)-1 + + paf2g = np.kron(np.eye(na), np.eye(pg+1, qa)).dot(paf2g) + pag2f = pag2f.dot(np.kron(np.eye(na), np.eye(qa, pg+1))) + pbf2g = np.kron(np.eye(nb), np.eye(pg+1, qb)).dot(pbf2g) + pbg2f = pbg2f.dot(np.kron(np.eye(nb), np.eye(qb, pg+1))) + + # Move to the glue space + xa = np.linspace(-1, 1, na + 1, endpoint=True) + xb = np.linspace(-1, 1, nb + 1, endpoint=True) + + # Glue-to-glue. + [vxi_g_dump, pa2g, pg2a, pb2g, pg2b] = make_projection_g2g_hr(pg, xa, xb) + + # Stack the operators. + pa2b = ((pbg2f.dot(pg2b)).dot(pa2g)).dot(paf2g) + pb2a = ((pag2f.dot(pg2a)).dot(pb2g)).dot(pbf2g) + + return pa2b, pb2a + + +def glue_to_dg(order, gg, dg_nodes): + + from modepy.quadrature.jacobi_gauss import legendre_gauss_lobatto_nodes + # Each element has order + 1 nodes. + nelem = int(dg_nodes.shape[0] / (order + 1)) + + # Create projection operators. + pglue2dg = np.zeros((dg_nodes.shape[0], (gg.shape[0]-1)*order)) + pdg2glue = np.zeros(((gg.shape[0]-1)*order, dg_nodes.shape[0])) + + for i_elem in range(0, nelem): + + a_elem = dg_nodes[i_elem*(order+1)] + b_elem = dg_nodes[i_elem*(order+1)+order] + + # Find the Legendre coefficients for coincident GG interval, using + # appropriate-order LGL quadrature. + quad_nodes = legendre_gauss_lobatto_nodes(order) + # Before transforming these, use them to get the weights. + quad_weights = lgl_weights(order+1, quad_nodes) + # Transform nodes to the appropriate interval. + for i in range(0, order + 1): + quad_nodes[i] = (quad_nodes[i] + 1) * (b_elem - a_elem) \ + / 2.0 + a_elem + # Transform the weights as well. + for i in range(0, order + 1): + quad_weights[i] = quad_weights[i] * (b_elem - a_elem) / 2.0 + # Verify that nodal points on the interface are already LGL points. + for i in range(0, order+1): + assert abs(quad_nodes[i] - dg_nodes[i_elem*(order+1)+i]) < 1e-12 + + # Do the quadrature. + # Coefficient loop. + for i in range(0, order): + # Quadrature point loop. + for k in range(0, order+1): + # Get Legendre polynomial of order i in + # interval, evaluated at quad point k. + poly_at_k = legendre_interval(i, a_elem, + b_elem, quad_nodes[k]) + + # Get Legendre coefficients. + i_op = i_elem*order + i + j_op = k + i_elem*(order+1) + pdg2glue[i_op][j_op] = ((2*i+1)/(b_elem - a_elem))\ + * quad_weights[k] * poly_at_k + + # Now that we have the projection from DG to glue, the projection + # from glue to DG comes from the H-compatibility condition + + # In our case, H is a repeated diagonal of the quad weights. + h = np.zeros((nelem*(order+1), nelem*(order+1))) + for i_elem in range(0, nelem): + for i in range(0, order+1): + h[i_elem*(order+1)+i][i_elem*(order+1)+i] = quad_weights[i] + + # Will need to invert this. + hinv = np.linalg.inv(h) + + # Now we need to also construct the per-interval mass matrix, + # taking advantage of the fact that the Legendre polynomials + # are orthogonal. + + m = np.zeros(((gg.shape[0]-1)*order, (gg.shape[0]-1)*order)) + m_small = np.zeros(order) + # Only same-order Legendre polynomials in the same interval + # will be nonzero, so the matrix is diagonal. + for i in range(0, order): + m_small[i] = (gg[1] - gg[0]) / (2.0 * i + 1) + for i in range(0, (gg.shape[0]-1)*order): + i_rem = i % order + m[i][i] = m_small[i_rem] + + # Finally, we can use compatibility to get the reverse projection. + pglue2dg = (m.dot(pdg2glue)).dot(hinv) + pglue2dg = pglue2dg.transpose() + + return pdg2glue, pglue2dg + + +def legendre_interval(order, a_int, b_int, point_k): + + # First check if the query point + # is within the interval. If it isn't, + # return 0. + if point_k >= a_int and point_k <= (b_int + 1e-12): + pass + else: + return 0 + + # With that out of the way, do some polynomial building + # using Bonnet's recursion formula, remembering to + # shift the Legendre polynomials to the interval + x_trans = (point_k - a_int) / (b_int - a_int) \ + - (b_int - point_k) / (b_int - a_int) + p_n_0 = 1 + p_n_1 = x_trans + + # If our input order is either of these, return them. + if order == 0: + return p_n_0 + elif order == 1: + return p_n_1 + + # Otherwise we enter recursion mode. + p_n_m1 = p_n_0 + p_n = p_n_1 + for i in range(1, order): + p_n_p1 = (1.0 / (i + 1)) * ((2*i + 1) * x_trans * p_n - i * p_n_m1) + p_n_m1 = p_n + p_n = p_n_p1 + + return p_n_p1 + + +def lgl_weights(n, nodes): + + weights = np.zeros(n) + for i in range(0, n): + weights[i] = (2.0 / (n*(n-1) + * (legendre_interval(n-1, -1, 1, nodes[i]))**2)) + + return weights diff --git a/test/q_2.txt b/test/q_2.txt new file mode 100644 index 000000000..e6ac0268e --- /dev/null +++ b/test/q_2.txt @@ -0,0 +1,6 @@ + 0.5000000000000000 + 0.1666666666666667 + 0.5000000000000000 + -0.1666666666666667 + 1.0000000000000000 + -0.3333333333333333 diff --git a/test/q_4.txt b/test/q_4.txt new file mode 100644 index 000000000..4b7cdf12b --- /dev/null +++ b/test/q_4.txt @@ -0,0 +1,96 @@ +-0.0416666666666667 +-0.0027777777777778 +0.0083333333333331 +0.0011904761904764 +0.5416666666666667 +0.1749999999999999 +-0.0083333333333331 +-0.0035714285714288 +0.5416666666666664 +-0.1749999999999997 +-0.0083333333333330 +0.0035714285714288 +-0.0416666666666666 +0.0027777777777777 +0.0083333333333332 +-0.0011904761904764 +1.2400202320584575 +0.0078507119289183 +0.0000000529778775 +-0.0000000352034607 +-0.0729337579017062 +0.0078478800177192 +-0.0000029202645105 +0.0000013131437493 +-0.0178405986993549 +0.0037727491934025 +0.0369615952094064 +0.0112488842016829 +0.0890096133619981 +0.2554389853632429 +0.0276743188599176 +-0.0165123473925149 +-0.2382554888193933 +0.3143058144807063 +-0.0201849049285142 +-0.0048899974275327 +0.3871547721536256 +-0.0850853038136043 +0.0000000347976171 +0.0000004882756223 +0.5704741035559282 +-0.0850856242077431 +0.0000013542349111 +-0.0000002055723806 +0.0683256772889949 +-0.0859279491160664 +-0.0050004014450767 +-0.0024166148237313 +-0.1013262031464566 +-0.0318424986928614 +-0.0075438526011718 +0.0031249334078735 +0.0753716501479071 +-0.1213255227420779 +-0.1027984715767137 +-0.0139918812156549 +0.2575512802532740 +0.0381320823634987 +-0.0000001583253617 +-0.0000012981662079 +0.1954343589159886 +0.0381363203630235 +-0.0000002527030111 +-0.0000009933207050 +0.3452893165945521 +0.0421803361278826 +-0.0208138136083886 +-0.0053811555335502 +0.1259767834841859 +-0.2124800017310330 +-0.0028189687190127 +0.0096801097150306 +0.0757482607519996 +-0.0398415082546706 +0.3060383673815450 +0.0441960895872428 +-0.1427996031213978 +0.0662630873587877 +0.0000000786593311 +0.0000005634968173 +0.1464943743090148 +0.0662607365784819 +-0.0000003956966040 +0.0000006636350413 +0.6413190735565835 +0.0678610437019278 +0.0032993990175575 +0.0025631851830881 +0.4703689973197340 +-0.0325261616434401 +0.0019559089340423 +-0.0041963299765795 +-0.1153828420639347 +0.0747250037529541 +-0.1296203385835838 +-0.0214066073755637 diff --git a/test/q_6.txt b/test/q_6.txt new file mode 100644 index 000000000..5ef09bd35 --- /dev/null +++ b/test/q_6.txt @@ -0,0 +1,324 @@ +0.0076388888888890 +0.0002976190476190 +-0.0015873015873015 +-0.0001322751322751 +0.0000330687830690 +0.0000030062530061 +-0.0645833333333338 +-0.0042658730158725 +0.0130952380952380 +0.0018518518518519 +-0.0000992063492062 +-0.0000150312650316 +0.5569444444444445 +0.1779761904761896 +-0.0115079365079365 +-0.0048941798941798 +0.0000661375661373 +0.0000300625300628 +0.5569444444444448 +-0.1779761904761902 +-0.0115079365079367 +0.0048941798941798 +0.0000661375661374 +-0.0000300625300626 +-0.0645833333333338 +0.0042658730158727 +0.0130952380952383 +-0.0018518518518517 +-0.0000992063492060 +0.0000150312650315 +0.0076388888888888 +-0.0002976190476190 +-0.0015873015873016 +0.0001322751322751 +0.0000330687830688 +-0.0000030062530062 +1.4268563455122170 +0.0382682224832046 +0.0062730242497421 +-0.0000007965781015 +-0.0000021799552952 +0.0000031243708247 +-0.1850971670076224 +0.0622469797290170 +0.0062744190375794 +0.0000002829327749 +-0.0000034557680218 +0.0000022076848745 +-0.0648981557548598 +0.0862272866341048 +0.0062769233641760 +0.0000035338545074 +0.0000014627798546 +-0.0000002921816912 +-0.1581346143995408 +0.1095251374450681 +0.0145814029712218 +0.0012769461829696 +-0.0005373905643547 +-0.0000766303009154 +-0.2594155138207218 +0.1426955295429118 +-0.0474531350117842 +-0.0155963333944991 +0.0006274045660508 +0.0002516855883443 +0.1957195334491137 +-0.2312509699825552 +-0.0353392218318904 +0.0181934735669320 +0.0005316092509378 +-0.0002603302379771 +-0.0095926586883630 +-0.1030241153958490 +0.0285461026469565 +-0.0024152925720598 +-0.0006067429675004 +0.0000840873716417 +0.0545622307097761 +-0.0567239766258600 +0.0293759307346910 +0.0019277778383208 +0.0002123757825960 +0.0000361878454639 +0.3392129585195110 +-0.1507565990718361 +0.0034987545283517 +0.0000000552211631 +0.0000013909237075 +-0.0000017845166932 +0.5611964896479856 +-0.1161082144029679 +0.0034983690188400 +-0.0000001858413508 +0.0000017049198556 +-0.0000006888382763 +0.1417146396683978 +-0.0814598481212354 +0.0034971885203090 +-0.0000010473718728 +-0.0000007041460643 +0.0000006574649649 +-0.0533525569340640 +-0.0467665678596720 +0.0022661137186922 +-0.0002451535334260 +0.0001236463156133 +0.0000205249914517 +0.0600295356643070 +-0.0127942396309124 +0.0120613485586357 +0.0030500892829083 +-0.0001574385845970 +-0.0000696747062982 +-0.0457592973586723 +0.0482087389842354 +0.0061406967811889 +-0.0041461246843002 +-0.0001141261666472 +0.0000726824264687 +0.0204208168834889 +0.0214306912296826 +-0.0065609219424412 +0.0006995593026804 +0.0002176232370328 +-0.0000179572367564 +-0.0234625860909541 +0.0298929334512873 +0.0262264569996840 +0.0037531278946722 +-0.0006161601985534 +-0.0000814417976086 +0.3748942611599946 +0.1668117925907492 +-0.0181769376775942 +0.0000016505102586 +-0.0000030182771242 +0.0000029438211654 +0.0180139527050320 +0.0857028815667877 +-0.0181784679818216 +0.0000009910467807 +-0.0000013440918030 +-0.0000013064971565 +0.0746238386423028 +0.0045891148305835 +-0.0181778229909092 +-0.0000026458877929 +0.0000006240076727 +-0.0000030215871136 +0.6650168973148122 +-0.0756868842396258 +-0.0284338587698686 +-0.0014908662915176 +0.0005533686188464 +0.0000648984500317 +0.0825405356562529 +-0.1679579805970867 +0.0464787505087289 +0.0180197770250175 +-0.0005991907450912 +-0.0002030388131742 +-0.1271337770068975 +0.2395425466352896 +0.0439065706339785 +-0.0189271275362692 +-0.0005931549885790 +0.0002081702648572 +-0.0554008864222944 +0.1345998216307505 +-0.0305376118824286 +0.0012238676357885 +0.0003578248053824 +-0.0000866448053757 +-0.0325548220492018 +0.0193656942897140 +-0.1241249974282292 +-0.0143897631796122 +0.0016891849092679 +0.0002284663733918 +-0.1024581261420299 +0.1257909677432488 +0.0015153781862069 +-0.0000008474008065 +0.0000007789434931 +-0.0000004244145945 +0.2304391611706006 +0.0891571891704643 +0.0015167407319010 +-0.0000009748827134 +-0.0000006240586523 +0.0000009923667029 +0.6440801174119887 +0.0525269565693532 +0.0015173223183447 +0.0000020585296105 +0.0000000627324804 +0.0000009084112755 +0.0820245802248009 +0.0158215065016135 +0.0054897695036996 +0.0007655707893505 +-0.0003292871572016 +-0.0000437729607812 +0.0981630408361174 +-0.0197771964741964 +-0.0240721736921666 +-0.0093092692657951 +0.0003701628210641 +0.0001405185367588 +0.0083452908862361 +-0.0989989738815951 +-0.0192029883069766 +0.0103760447110098 +0.0003403703897566 +-0.0001440922422889 +-0.0288188490304330 +-0.0859982980988029 +0.0177618414433881 +-0.0004103013388925 +-0.0004612200444481 +0.0000380919055893 +0.0682247846427184 +-0.0859464570522154 +-0.0628137261290580 +-0.0083111162477239 +0.0010391751207186 +0.0001087956757791 +-0.1011503113778986 +-0.0162476416019815 +0.0011983954402854 +-0.0000000588179514 +-0.0000005947179354 +0.0000005482563072 +0.0719904403159826 +-0.0019796478930281 +0.0011977824458625 +0.0000012499070772 +0.0000000613542417 +0.0000003740606944 +0.0832538157106969 +0.0122861500314558 +0.0011977298182762 +0.0000001277266936 +0.0000003015490923 +0.0000002630332989 +0.6566754683652896 +0.0262616468261796 +0.0029451944128836 +0.0001330300433277 +-0.0000334377442731 +-0.0000040227100399 +0.1400781774850962 +0.0446073035815920 +-0.0109409719384444 +-0.0017177157691836 +0.0000669932825490 +0.0000159076164692 +0.1257481536103612 +-0.1146025911016925 +-0.0037530097519190 +0.0028742337159273 +0.0000373772981286 +-0.0000201476048079 +0.0932564138521604 +0.0035389065504577 +0.0035428896959126 +-0.0013480068067448 +0.0002762792753176 +0.0000358035921204 +-0.0698521579616875 +0.1098332091938433 +0.2163649169891393 +0.0260966655396425 +-0.0029253556960296 +-0.0003337393398058 +0.0607557809217103 +-0.0477299636863383 +0.0015660752104962 +0.0000002406122343 +0.0000002217587330 +-0.0000003222125190 +-0.0834412984870233 +-0.0205232662007936 +0.0015660005245034 +-0.0000003778473675 +0.0000002791480044 +-0.0000004851896970 +-0.0969053718693941 +0.0066836149296925 +0.0015657755275141 +-0.0000006607220888 +-0.0000002243360089 +-0.0000002883209355 +-0.0012435327466196 +0.0335825687247331 +0.0021479970891172 +-0.0000648055397491 +0.0000557310904477 +0.0000097726914592 +0.7438158486829837 +0.0650351693707065 +-0.0025803264052022 +0.0007626503461785 +-0.0000566482134126 +-0.0000315812989750 +0.4503079338353286 +-0.0897576916820822 +0.0022852763449330 +-0.0004173965860673 +-0.0000920770941274 +0.0000292294123977 +-0.0954825204170535 +0.0253564541898826 +0.0054354864750859 +-0.0009460553470582 +-0.0000801889312214 +-0.0000148890576181 +0.0221931600800686 +-0.0292017286479252 +-0.0875580840377385 +-0.0100077065308661 +0.0011248933398002 +0.0001230164560477 diff --git a/test/sbp_operators.py b/test/sbp_operators.py new file mode 100644 index 000000000..fdcbbe089 --- /dev/null +++ b/test/sbp_operators.py @@ -0,0 +1,163 @@ +import numpy as np + + +# Standard FD method: second-order +def sbp21(n): + p = np.zeros((n, n, )) + q = np.zeros((n, n, )) + + # norm matrix + p[0, 0] = 0.5 + p[n-1, n-1] = 0.5 + for i in range(1, n-1): + p[i, i] = 1.0 + + # now the q matrix + q[0, 1] = 0.5 + q[n-1, n-2] = -0.5 + for i in range(1, n-1): + q[i, i-1] = -0.5 + q[i, i+1] = 0.5 + + return p, q + + +# Standard FD method: third-order +def sbp42(n): + p = np.zeros((n, n, )) + q = np.zeros((n, n, )) + + # norm matrix + # upper subblock + p[0, 0] = 17.0/48.0 + p[1, 1] = 59.0/48.0 + p[2, 2] = 43.0/48.0 + p[3, 3] = 49.0/48.0 + # lower subblock + p[n-1, n-1] = 17.0/48.0 + p[n-2, n-2] = 59.0/48.0 + p[n-3, n-3] = 43.0/48.0 + p[n-4, n-4] = 49.0/48.0 + for i in range(4, n-4): + p[i, i] = 1.0 + + # now the q matrix + # upper subblock + q[0, 1] = 59.0/96.0 + q[0, 2] = -1.0/12.0 + q[0, 3] = -1.0/32.0 + q[1, 2] = 59.0/96.0 + q[2, 3] = 59.0/96.0 + q[1, 0] = -q[0, 1] + q[2, 0] = -q[0, 2] + q[2, 1] = -q[1, 2] + q[2, 4] = -1.0/12.0 + q[3, 0] = -q[0, 3] + q[3, 2] = -q[2, 3] + q[3, 4] = 2.0/3.0 + q[3, 5] = -1.0/12.0 + # lower subblock + q[n-1, n-2] = -q[0, 1] + q[n-1, n-3] = -q[0, 2] + q[n-1, n-4] = -q[0, 3] + q[n-2, n-1] = -q[1, 0] + q[n-2, n-3] = -q[1, 2] + q[n-3, n-1] = -q[2, 0] + q[n-3, n-2] = -q[2, 1] + q[n-3, n-4] = -q[2, 3] + q[n-3, n-5] = 1.0/12.0 + q[n-4, n-1] = -q[3, 0] + q[n-4, n-3] = -q[3, 2] + q[n-4, n-5] = -2.0/3.0 + q[n-4, n-6] = 1.0/12.0 + for i in range(4, n-4): + q[i, i-2] = 1.0/12.0 + q[i, i-1] = -2.0/3.0 + q[i, i+1] = 2.0/3.0 + q[i, i+2] = -1.0/12.0 + + return p, q + + +# Standard FD method: fourth-order +def sbp63(n): + p = np.zeros((n, n, )) + q = np.zeros((n, n, )) + + # norm matrix + # upper subblock + p[0, 0] = 13649.0/43200.0 + p[1, 1] = 12013.0/8640.0 + p[2, 2] = 2711.0/4320.0 + p[3, 3] = 5359.0/4320.0 + p[4, 4] = 7877.0/8640.0 + p[5, 5] = 43801.0/43200.0 + # lower subblock + p[n-1, n-1] = 13649.0/43200.0 + p[n-2, n-2] = 12013.0/8640.0 + p[n-3, n-3] = 2711.0/4320.0 + p[n-4, n-4] = 5359.0/4320.0 + p[n-5, n-5] = 7877.0/8640.0 + p[n-6, n-6] = 43801.0/43200.0 + for i in range(6, n-6): + p[i, i] = 1.0 + + # now the q matrix + # upper subblock + q[0, 1] = 127/211 + q[0, 2] = 35/298 + q[0, 3] = -32/83 + q[0, 4] = 89/456 + q[0, 5] = -2/69 + + q[1, 0] = -127/211 + q[1, 2] = -1/167 + q[1, 3] = 391/334 + q[1, 4] = -50/71 + q[1, 5] = 14/99 + + q[2, 0] = -35/298 + q[2, 1] = 1/167 + q[2, 3] = -34/79 + q[2, 4] = 90/113 + q[2, 5] = -69/271 + + q[3, 0] = 32/83 + q[3, 1] = -391/334 + q[3, 2] = 34/79 + q[3, 4] = 6/25 + q[3, 5] = 31/316 + q[3, 6] = 1/60 + + q[4, 0] = -89/456 + q[4, 1] = 50/71 + q[4, 2] = -90/113 + q[4, 3] = -6/25 + q[4, 5] = 37/56 + q[4, 6] = -3/20 + q[4, 7] = 1/60 + + q[5, 0] = 2/69 + q[5, 1] = -14/99 + q[5, 2] = 69/271 + q[5, 3] = -31/316 + q[5, 4] = -37/56 + q[5, 6] = 3/4 + q[5, 7] = -3/20 + q[5, 8] = 1/60 + + # lower subblock + for i in range(0, 6): + for j in range(0, 9): + q[n-1-i, n-1-j] = -q[i, j] + + # interior + for i in range(6, n-6): + q[i, i-3] = -1.0/60.0 + q[i, i-2] = 3.0/20.0 + q[i, i-1] = -3.0/4.0 + q[i, i+1] = 3.0/4.0 + q[i, i+2] = -3.0/20.0 + q[i, i+3] = 1.0/60.0 + + return p, q diff --git a/test/test_grudge.py b/test/test_grudge.py index f6e4ad9b8..e72e30fa2 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -819,6 +819,374 @@ def rhs(t, u): else: assert eoc_rec.order_estimate() > order + +@pytest.mark.parametrize("order", [4]) +@pytest.mark.parametrize("order_sbp", [6]) +@pytest.mark.parametrize("spacing_factor", [2]) +@pytest.mark.parametrize("c", [np.array([0.27, 0.31]), np.array([-0.27, 0.31])]) +def test_convergence_sbp_advec(actx_factory, order, order_sbp, spacing_factor, c, + visualize=False): + + actx = actx_factory() + + # Domain: x = [-1,1], y = [-1,1] + # SBP Subdomain: x = [-1,0], y = [0,1] (structured mesh) + # DG Subdomain: x = [0,1], y = [0,1] (unstructured mesh) + + # DG Half. + dim = 2 + + from pytools.convergence import EOCRecorder + eoc_rec = EOCRecorder() + from grudge import sym + from sbp_operators import (sbp21, sbp42, sbp63) + from projection import sbp_dg_projection + + nelems = [20, 40, 60] + for nelem in nelems: + + mesh = mgen.generate_regular_rect_mesh(a=(0, -1), b=(1, 1), + n=(nelem, nelem), order=order, + boundary_tag_to_face={ + "btag_sbp": ["-x"], + "btag_std": ["+y", "+x", + "-y"]}) + + # Check to make sure this actually worked. + from meshmode.mesh import is_boundary_tag_empty + assert not is_boundary_tag_empty(mesh, "btag_sbp") + assert not is_boundary_tag_empty(mesh, "btag_std") + + # Check this new isolating discretization, as well as the inverse. + dcoll = DiscretizationCollection(actx, mesh, order=order) + sbp_bdry_discr = dcoll.discr_from_dd(dof_desc.DTAG_BOUNDARY("btag_sbp")) + sbp_bdry_discr.cl_context = actx + + dt_factor = 4 + h = 1/nelem + + norm_c = la.norm(c) + + flux_type = "upwind" + + def f(x): + return actx.np.sin(10*x) + + def u_analytic(x, t=0): + return f(-c.dot(x)/norm_c+t*norm_c) + + from grudge.models.advection import WeakAdvectionSBPOperator + + std_tag = dof_desc.DTAG_BOUNDARY("btag_std") + adv_operator = WeakAdvectionSBPOperator(dcoll, c, + inflow_u=lambda t: u_analytic( + thaw(dcoll.nodes(dd=std_tag), + actx), + t=t + ), + flux_type=flux_type) + + nodes = thaw(dcoll.nodes(), actx) + u = u_analytic(nodes, t=0) + + # Count the number of DG nodes - will need this later + ngroups = len(mesh.groups) + nnodes_grp = np.ones(ngroups) + for i in range(0, dim): + for j in range(0, ngroups): + nnodes_grp[j] = nnodes_grp[j] * \ + mesh.groups[j].nodes.shape[i*-1 - 1] + + nnodes = int(sum(nnodes_grp)) + + final_time = 0.2 + + dt = dt_factor * h/order**2 + nsteps = (final_time // dt) + 1 + dt = final_time/nsteps + 1e-15 + + # SBP Half. + + # First, need to set up the structured mesh. + n_sbp_x = nelem + n_sbp_y = nelem*spacing_factor + x_sbp = np.linspace(-1, 0, n_sbp_x, endpoint=True) + y_sbp = np.linspace(-1, 1, n_sbp_y, endpoint=True) + dx = x_sbp[1] - x_sbp[0] + dy = y_sbp[1] - y_sbp[0] + + # Set up solution vector: + # For now, timestep is the same as DG. + u_sbp = actx.np.zeros(int(n_sbp_x*n_sbp_y)) + + # Initial condition + for j in range(0, n_sbp_y): + for i in range(0, n_sbp_x): + u_sbp[i + j*(n_sbp_x)] = np.sin(10*(-c.dot( + [x_sbp[i], + y_sbp[j]])/norm_c)) + + # obtain P and Q + if order_sbp == 2: + [p_x, q_x] = sbp21(n_sbp_x) + [p_y, q_y] = sbp21(n_sbp_y) + elif order_sbp == 4: + [p_x, q_x] = sbp42(n_sbp_x) + [p_y, q_y] = sbp42(n_sbp_y) + elif order_sbp == 6: + [p_x, q_x] = sbp63(n_sbp_x) + [p_y, q_y] = sbp63(n_sbp_y) + + tau_l = 1 + tau_r = 1 + + # for the boundaries + el_x = np.zeros(n_sbp_x) + er_x = np.zeros(n_sbp_x) + el_x[0] = 1 + er_x[n_sbp_x-1] = 1 + e_l_matx = np.zeros((n_sbp_x, n_sbp_x,)) + e_r_matx = np.zeros((n_sbp_x, n_sbp_x,)) + + for i in range(0, n_sbp_x): + for j in range(0, n_sbp_x): + e_l_matx[i, j] = el_x[i]*el_x[j] + e_r_matx[i, j] = er_x[i]*er_x[j] + + el_y = np.zeros(n_sbp_y) + er_y = np.zeros(n_sbp_y) + el_y[0] = 1 + er_y[n_sbp_y-1] = 1 + e_l_maty = np.zeros((n_sbp_y, n_sbp_y,)) + e_r_maty = np.zeros((n_sbp_y, n_sbp_y,)) + + for i in range(0, n_sbp_y): + for j in range(0, n_sbp_y): + e_l_maty[i, j] = el_y[i]*el_y[j] + e_r_maty[i, j] = er_y[i]*er_y[j] + + # construct the spatial operators + d_x = np.linalg.inv(dx*p_x).dot(q_x - 0.5*e_l_matx + 0.5*e_r_matx) + d_y = np.linalg.inv(dy*p_y).dot(q_y - 0.5*e_l_maty + 0.5*e_r_maty) + + # for the boundaries + c_l_x = np.kron(tau_l, (np.linalg.inv(dx*p_x).dot(el_x))) + c_r_x = np.kron(tau_r, (np.linalg.inv(dx*p_x).dot(er_x))) + c_l_y = np.kron(tau_l, (np.linalg.inv(dy*p_y).dot(el_y))) + c_r_y = np.kron(tau_r, (np.linalg.inv(dy*p_y).dot(er_y))) + + # For speed... + dudx_mat = -np.kron(np.eye(n_sbp_y), d_x) + dudy_mat = -np.kron(d_y, np.eye(n_sbp_x)) + + # Number of nodes in our SBP-DG boundary discretization + sbp_nodes_y = thaw(sbp_bdry_discr.nodes(), actx)[1] + # When projecting, we use nodes sorted in y, but we will have to unsort + # afterwards to make sure projected solution is injected into DG BC + # in the correct way. + nodesort = np.argsort(sbp_nodes_y.get()) + nodesortlist = nodesort.tolist() + rangex = np.array(range(sbp_nodes_y.get().shape[0])) + unsort_args = [nodesortlist.index(x) for x in rangex] + + west_nodes = np.sort(np.array(sbp_nodes_y.get())) + + # Make element-aligned glue grid. + dg_side_gg = np.zeros(int(west_nodes.shape[0]/(order+1))+1) + counter = 0 + for i in range(0, west_nodes.shape[0]): + if i % (order+1) == 0: + dg_side_gg[counter] = west_nodes[i] + counter += 1 + + dg_side_gg[-1] = west_nodes[-1] + n_west_elements = int(west_nodes.shape[0] / (order + 1)) + sbp2dg, dg2sbp = sbp_dg_projection(n_sbp_y-1, n_west_elements, + order_sbp, order, dg_side_gg, + west_nodes) + + # Get mapping for western face + base_nodes = thaw(dcoll._volume_discr.nodes()) + nsbp_nodes = sbp_bdry_discr.nnodes + nodes_per_element = mesh.groups[0].nodes.shape[2] + west_indices = np.zeros(nsbp_nodes) + count = 0 + # Sweep through first block of indices in the box. + for i in range(0, (nelem-1)*2*nodes_per_element): + if base_nodes[0][i] < 1e-10: + # Make sure we're actually at the edge faces. + if i % (2*nodes_per_element) < nodes_per_element: + west_indices[count] = i + count += 1 + + def rhs(t, u): + # Initialize the entire RHS to 0. + rhs_out = actx.np.zeros(int(n_sbp_x*n_sbp_y) + int(nnodes)) + + # Fill the first part with the SBP half of the domain. + + # Pull the SBP vector out of device array for now. + u_sbp_ts = u[0:int(n_sbp_x*n_sbp_y)].get() + + dudx = np.zeros((n_sbp_x*n_sbp_y)) + dudy = np.zeros((n_sbp_x*n_sbp_y)) + + dudx = dudx_mat.dot(u_sbp_ts) + dudy = dudy_mat.dot(u_sbp_ts) + + # Boundary condition + dl_x = np.zeros(n_sbp_x*n_sbp_y) + dr_x = np.zeros(n_sbp_x*n_sbp_y) + dl_y = np.zeros(n_sbp_x*n_sbp_y) + dr_y = np.zeros(n_sbp_x*n_sbp_y) + + # Pull DG solution at western face to project. + u_dg_ts = u[int(n_sbp_x*n_sbp_y):].get() + + dg_west = np.zeros(nsbp_nodes) + for i in range(0, nsbp_nodes): + dg_west[i] = u_dg_ts[int(west_indices[i])] + + # Project DG to SBP: + dg_proj = dg2sbp.dot(dg_west) + + # Need to fill this by looping through each segment. + # X-boundary conditions: + for j in range(0, n_sbp_y): + u_bcx = u_sbp_ts[j*n_sbp_x:((j+1)*n_sbp_x)] + v_l_x = np.transpose(el_x).dot(u_bcx) + v_r_x = np.transpose(er_x).dot(u_bcx) + left_bcx = np.sin(10*(-c.dot( + [x_sbp[0], y_sbp[j]])/norm_c + + norm_c*t)) + # Incorporate DG here. + right_bcx = dg_proj[j] + dl_xbc = c_l_x*(v_l_x - left_bcx) + dr_xbc = c_r_x*(v_r_x - right_bcx) + dl_x[j*n_sbp_x:((j+1)*n_sbp_x)] = dl_xbc + dr_x[j*n_sbp_x:((j+1)*n_sbp_x)] = dr_xbc + # Y-boundary conditions: + for i in range(0, n_sbp_x): + u_bcy = u_sbp_ts[i::n_sbp_x] + v_l_y = np.transpose(el_y).dot(u_bcy) + v_r_y = np.transpose(er_y).dot(u_bcy) + left_bcy = np.sin(10*(-c.dot( + [x_sbp[i], y_sbp[0]])/norm_c + norm_c*t)) + right_bcy = np.sin(10*(-c.dot( + [x_sbp[i], + y_sbp[n_sbp_y-1]])/norm_c + norm_c*t)) + dl_ybc = c_l_y*(v_l_y - left_bcy) + dr_ybc = c_r_y*(v_r_y - right_bcy) + dl_y[i::n_sbp_x] = dl_ybc + dr_y[i::n_sbp_x] = dr_ybc + + # Add these at each point on the SBP half to get the SBP RHS. + rhs_sbp = c[0]*dudx + c[1]*dudy - dl_x - dr_x - dl_y - dr_y + + # Now pop this back into the device RHS vector. + rhs_sbp_dev = actx.np.zeros((n_sbp_x*n_sbp_y,)) + rhs_sbp_dev = rhs_sbp + + rhs_out[0:int(n_sbp_x*n_sbp_y)] = rhs_sbp_dev + + sbp_east = np.zeros(n_sbp_y) + # Pull SBP domain values off of east face. + counter = 0 + for i in range(0, n_sbp_x*n_sbp_y): + if i == n_sbp_x - 1: + sbp_east[counter] = u_sbp_ts[i] + counter += 1 + elif i % n_sbp_x == n_sbp_x - 1: + sbp_east[counter] = u_sbp_ts[i] + counter += 1 + + # Projection from SBP to DG is now a two-step process. + # First: SBP-to-DG. + sbp_proj = sbp2dg.dot(sbp_east) + # Second: Fix the ordering. + sbp_proj = sbp_proj[unsort_args] + + # Grudge DG RHS. + # Critical step - now need to apply projected SBP state to the + # proper nodal locations in u_dg. + rhs_out[int(n_sbp_x*n_sbp_y):] = adv_operator.operator( + t=t, + u=u[int(n_sbp_x*n_sbp_y):], + sbp_tag=sym.DTAG_BOUNDARY("btag_sbp"), std_tag=std_tag, + state_from_sbp=sbp_proj) + + return rhs_out + + # Timestepper. + from grudge.shortcuts import set_up_rk4 + + # Make a combined u with the SBP and the DG parts. + u_comb = actx.np.zeros(int(n_sbp_x*n_sbp_y) + nnodes) + u_comb[0:int(n_sbp_x*n_sbp_y)] = u_sbp + for i in range(int(n_sbp_x*n_sbp_y), int(n_sbp_x*n_sbp_y) + nnodes): + u_comb[i] = u[i - int(n_sbp_x*n_sbp_y)] + dt_stepper = set_up_rk4("u", dt, u_comb, rhs) + + from grudge.shortcuts import make_visualizer + vis = make_visualizer(dcoll, vis_order=order) + + step = 0 + + # Create mesh for structured grid output + sbp_mesh = np.zeros((2, n_sbp_y, n_sbp_x,)) + for j in range(0, n_sbp_y): + sbp_mesh[0, j, :] = x_sbp + for i in range(0, n_sbp_x): + sbp_mesh[1, :, i] = y_sbp + + for event in dt_stepper.run(t_end=final_time): + if isinstance(event, dt_stepper.StateComputed): + + step += 1 + + last_t = event.t + u_sbp = event.state_component[0:int(n_sbp_x*n_sbp_y)] + u_dg = event.state_component[int(n_sbp_x*n_sbp_y):] + + error_l2_dg = op.norm( + dcoll, + u_dg - u_analytic(nodes, t=last_t), + 2 + ) + + sbp_error = np.zeros((n_sbp_x*n_sbp_y)) + error_l2_sbp = 0 + for j in range(0, n_sbp_y): + for i in range(0, n_sbp_x): + sbp_error[i + j*n_sbp_x] = \ + u_sbp[i + j*n_sbp_x].get() - \ + np.sin(10*(-c.dot([x_sbp[i], y_sbp[j]]) + / norm_c + last_t*norm_c)) + error_l2_sbp = error_l2_sbp + \ + dx*dy*(sbp_error[i + j*n_sbp_x]) ** 2 + + error_l2_sbp = np.sqrt(error_l2_sbp) + + # Write out the DG data + if visualize: + vis.write_vtk_file("eoc_dg-%s-%04d.vtu" % + (nelem, step), + [("u", u_dg)], overwrite=True) + + # Write out the SBP data. + from pyvisfile.vtk import write_structured_grid + if visualize: + filename = "eoc_sbp_%s-%04d.vts" % (nelem, step) + write_structured_grid(filename, sbp_mesh, + point_data=[("u", u_sbp.get())]) + + if c[0] > 0: + eoc_rec.add_data_point(h, error_l2_dg) + else: + eoc_rec.add_data_point(h, error_l2_sbp) + + assert eoc_rec.order_estimate() > (order_sbp/2 + 1) * 0.95 + # }}} diff --git a/test/test_sbp_dg.py b/test/test_sbp_dg.py new file mode 100644 index 000000000..09eec6245 --- /dev/null +++ b/test/test_sbp_dg.py @@ -0,0 +1,393 @@ +# Copyright (C) 2020 Cory Mikida +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU 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 General Public License +# along with this program. If not, see . + + +import numpy as np +import numpy.linalg as la + +from grudge.array_context import PytestPyOpenCLArrayContextFactory +from arraycontext import pytest_generate_tests_for_array_contexts +pytest_generate_tests = pytest_generate_tests_for_array_contexts( + [PytestPyOpenCLArrayContextFactory]) + +from arraycontext.container.traversal import thaw +import meshmode.mesh.generation as mgen + +from grudge import DiscretizationCollection + +import grudge.op as op +import grudge.dof_desc as dof_desc +from sbp_operators import (sbp21, sbp42, sbp63) +from projection import sbp_dg_projection + + +import pytest + +import logging + +logger = logging.getLogger(__name__) + + +# Domain: x = [-1,1], y = [-1,1] +# SBP Subdomain: x = [-1,0], y = [0,1] (structured mesh) +# DG Subdomain: x = [0,1], y = [0,1] (unstructured mesh) + +def test_sbp_dg(actx_factory, write_output=True, order=4): + actx = actx_factory() + + # DG Half. + dim = 2 + + from grudge import sym + + nelem_x = 20 + nelem_y = 20 + mesh = mgen.generate_regular_rect_mesh(a=(0, -1), b=(1, 1), + n=(nelem_x, nelem_y), order=order, + boundary_tag_to_face={"btag_sbp": ["-x"], + "btag_std": ["+y", "+x", "-y"]}) + + # Check to make sure this actually worked. + from meshmode.mesh import is_boundary_tag_empty + assert not is_boundary_tag_empty(mesh, "btag_sbp") + assert not is_boundary_tag_empty(mesh, "btag_std") + + # Check this new isolating discretization, as well as the inverse. + dcoll = DiscretizationCollection(actx, mesh, order=order) + sbp_bdry_discr = dcoll.discr_from_dd(dof_desc.DTAG_BOUNDARY("btag_sbp")) + sbp_bdry_discr.cl_context = actx + + # Visualize our new boundary discretization + # from grudge.shortcuts import make_visualizer + # vis = make_visualizer(sbp_bdry_discr, vis_order=order) + # vis.write_vtk_file("bdry_disc.vtu", + # [], overwrite=True) + + std_bdry_discr = dcoll.discr_from_dd(dof_desc.DTAG_BOUNDARY("btag_std")) + std_bdry_discr.cl_context = actx + + # Visualize our new boundary discretization + # vis = make_visualizer(std_bdry_discr, vis_order=order) + # vis.write_vtk_file("std_bdry_disc.vtu", + # [], overwrite=True) + + dt_factor = 4 + h = 1/40 + + c = np.array([0.1, 0.1]) + norm_c = la.norm(c) + + flux_type = "upwind" + + def f(x): + return actx.np.sin(10*x) + + def u_analytic(x, t=0): + return f(-c.dot(x)/norm_c+t*norm_c) + + from grudge.models.advection import WeakAdvectionSBPOperator + + std_tag = dof_desc.DTAG_BOUNDARY("btag_std") + from meshmode.mesh import BTAG_ALL + + adv_operator = WeakAdvectionSBPOperator(dcoll, c, + inflow_u=lambda t: u_analytic( + thaw(dcoll.nodes(dd=std_tag), + actx), + t=t + ), + flux_type=flux_type) + + nodes = thaw(dcoll.nodes(), actx) + u = u_analytic(nodes, t=0) + + # Count the number of DG nodes - will need this later + ngroups = len(mesh.groups) + nnodes_grp = np.ones(ngroups) + for i in range(0, dim): + for j in range(0, ngroups): + nnodes_grp[j] = nnodes_grp[j]*mesh.groups[j].nodes.shape[i*-1 - 1] + + nnodes = int(sum(nnodes_grp)) + + final_time = 0.2 + + dt = dt_factor * h/order**2 + nsteps = (final_time // dt) + 1 + dt = final_time/nsteps + 1e-15 + + # SBP Half. + + # First, need to set up the structured mesh. + n_sbp_x = 20 + n_sbp_y = 40 + x_sbp = np.linspace(-1, 0, n_sbp_x, endpoint=True) + y_sbp = np.linspace(-1, 1, n_sbp_y, endpoint=True) + dx = x_sbp[1] - x_sbp[0] + dy = y_sbp[1] - y_sbp[0] + + # Set up solution vector: + # For now, timestep is the same as DG. + u_sbp = np.zeros(int(n_sbp_x*n_sbp_y)) + + # Initial condition + for j in range(0, n_sbp_y): + for i in range(0, n_sbp_x): + u_sbp[i + j*(n_sbp_x)] = np.sin(10*(-c.dot( + [x_sbp[i], y_sbp[j]])/norm_c)) + + # obtain P and Q + order_sbp = 4 + if order_sbp == 2: + [p_x, q_x] = sbp21(n_sbp_x) + [p_y, q_y] = sbp21(n_sbp_y) + elif order_sbp == 4: + [p_x, q_x] = sbp42(n_sbp_x) + [p_y, q_y] = sbp42(n_sbp_y) + elif order_sbp == 6: + [p_x, q_x] = sbp63(n_sbp_x) + [p_y, q_y] = sbp63(n_sbp_y) + + tau_l = 1 + tau_r = 1 + + # for the boundaries + el_x = np.zeros(n_sbp_x) + er_x = np.zeros(n_sbp_x) + el_x[0] = 1 + er_x[n_sbp_x-1] = 1 + e_l_matx = np.zeros((n_sbp_x, n_sbp_x,)) + e_r_matx = np.zeros((n_sbp_x, n_sbp_x,)) + + for i in range(0, n_sbp_x): + for j in range(0, n_sbp_x): + e_l_matx[i, j] = el_x[i]*el_x[j] + e_r_matx[i, j] = er_x[i]*er_x[j] + + el_y = np.zeros(n_sbp_y) + er_y = np.zeros(n_sbp_y) + el_y[0] = 1 + er_y[n_sbp_y-1] = 1 + e_l_maty = np.zeros((n_sbp_y, n_sbp_y,)) + e_r_maty = np.zeros((n_sbp_y, n_sbp_y,)) + + for i in range(0, n_sbp_y): + for j in range(0, n_sbp_y): + e_l_maty[i, j] = el_y[i]*el_y[j] + e_r_maty[i, j] = er_y[i]*er_y[j] + + # construct the spatial operators + d_x = np.linalg.inv(dx*p_x).dot(q_x - 0.5*e_l_matx + 0.5*e_r_matx) + d_y = np.linalg.inv(dy*p_y).dot(q_y - 0.5*e_l_maty + 0.5*e_r_maty) + + # for the boundaries + c_l_x = np.kron(tau_l, (np.linalg.inv(dx*p_x).dot(el_x))) + c_r_x = np.kron(tau_r, (np.linalg.inv(dx*p_x).dot(er_x))) + c_l_y = np.kron(tau_l, (np.linalg.inv(dy*p_y).dot(el_y))) + c_r_y = np.kron(tau_r, (np.linalg.inv(dy*p_y).dot(er_y))) + + # For speed... + dudx_mat = -np.kron(np.eye(n_sbp_y), d_x) + dudy_mat = -np.kron(d_y, np.eye(n_sbp_x)) + + # Number of nodes in our SBP-DG boundary discretization + from meshmode.dof_array import flatten + sbp_nodes_y = flatten(thaw(sbp_bdry_discr.nodes(), actx)[1]) + # When projecting, we use nodes sorted in y, but we will have to unsort + # afterwards to make sure projected solution is injected into DG BC + # in the correct way. + nodesort = np.argsort(sbp_nodes_y) + nodesortlist = nodesort.tolist() + rangex = np.array(range(sbp_nodes_y.shape[0])) + unsort_args = [nodesortlist.index(x) for x in rangex] + + west_nodes = np.sort(np.array(sbp_nodes_y)) + + # Make element-aligned glue grid. + dg_side_gg = np.zeros(int(west_nodes.shape[0]/(order+1))+1) + counter = 0 + for i in range(0, west_nodes.shape[0]): + west_nodes[i] = west_nodes[i].get() + if i % (order+1) == 0: + dg_side_gg[counter] = west_nodes[i] + counter += 1 + + dg_side_gg[-1] = west_nodes[-1] + n_west_elements = int(west_nodes.shape[0] / (order + 1)) + sbp2dg, dg2sbp = sbp_dg_projection(n_sbp_y-1, n_west_elements, order_sbp, + order, dg_side_gg, west_nodes) + + def rhs(t, u): + # Initialize the entire RHS to 0. + rhs_out = np.zeros(int(n_sbp_x*n_sbp_y) + int(nnodes)) + + # Fill the first part with the SBP half of the domain. + + # Pull the SBP vector out of device array for now. + u_sbp_ts = u[0:int(n_sbp_x*n_sbp_y)] + + dudx = np.zeros((n_sbp_x*n_sbp_y)) + dudy = np.zeros((n_sbp_x*n_sbp_y)) + + dudx = dudx_mat.dot(u_sbp_ts) + dudy = dudy_mat.dot(u_sbp_ts) + + # Boundary condition + dl_x = np.zeros(n_sbp_x*n_sbp_y) + dr_x = np.zeros(n_sbp_x*n_sbp_y) + dl_y = np.zeros(n_sbp_x*n_sbp_y) + dr_y = np.zeros(n_sbp_x*n_sbp_y) + + # Need to fill this by looping through each segment. + # X-boundary conditions: + for j in range(0, n_sbp_y): + u_bcx = u_sbp_ts[j*n_sbp_x:((j+1)*n_sbp_x)] + v_l_x = np.transpose(el_x).dot(u_bcx) + v_r_x = np.transpose(er_x).dot(u_bcx) + left_bcx = np.sin(10*(-c.dot( + [x_sbp[0], y_sbp[j]])/norm_c + norm_c*t)) + right_bcx = np.sin(10*(-c.dot( + [x_sbp[n_sbp_x-1], + y_sbp[j]])/norm_c + norm_c*t)) + dl_xbc = c_l_x*(v_l_x - left_bcx) + dr_xbc = c_r_x*(v_r_x - right_bcx) + dl_x[j*n_sbp_x:((j+1)*n_sbp_x)] = dl_xbc + dr_x[j*n_sbp_x:((j+1)*n_sbp_x)] = dr_xbc + # Y-boundary conditions: + for i in range(0, n_sbp_x): + u_bcy = u_sbp_ts[i::n_sbp_x] + v_l_y = np.transpose(el_y).dot(u_bcy) + v_r_y = np.transpose(er_y).dot(u_bcy) + left_bcy = np.sin(10*(-c.dot( + [x_sbp[i], y_sbp[0]])/norm_c + norm_c*t)) + right_bcy = np.sin(10*(-c.dot( + [x_sbp[i], + y_sbp[n_sbp_y-1]])/norm_c + norm_c*t)) + dl_ybc = c_l_y*(v_l_y - left_bcy) + dr_ybc = c_r_y*(v_r_y - right_bcy) + dl_y[i::n_sbp_x] = dl_ybc + dr_y[i::n_sbp_x] = dr_ybc + + # Add these at each point on the SBP half to get the SBP RHS. + rhs_sbp = c[0]*dudx + c[1]*dudy - dl_x - dr_x - dl_y - dr_y + + # Now pop this back into the device RHS vector. + # rhs_sbp_dev = actx.np.zeros((n_sbp_x*n_sbp_y,)) + rhs_sbp_dev = rhs_sbp + + rhs_out[0:int(n_sbp_x*n_sbp_y)] = rhs_sbp_dev + + sbp_east = np.zeros(n_sbp_y) + # Pull SBP domain values off of east face. + counter = 0 + for i in range(0, n_sbp_x*n_sbp_y): + if i == n_sbp_x - 1: + sbp_east[counter] = u_sbp_ts[i] + counter += 1 + elif i % n_sbp_x == n_sbp_x - 1: + sbp_east[counter] = u_sbp_ts[i] + counter += 1 + + # Projection from SBP to DG is now a two-step process. + # First: SBP-to-DG. + sbp_proj = sbp2dg.dot(sbp_east) + # Second: Fix the ordering. + sbp_proj = sbp_proj[unsort_args] + + # Grudge DG RHS. + # Critical step - now need to apply projected SBP state to the + # proper nodal locations in u_dg. + rhs_out[int(n_sbp_x*n_sbp_y):] = adv_operator.operator( + t=t, + u=u[int(n_sbp_x*n_sbp_y):], + state_from_sbp=sbp_proj, + sbp_tag=dof_desc.DTAG_BOUNDARY("btag_sbp"), std_tag=std_tag) + + return rhs_out + + # Timestepper. + from grudge.shortcuts import set_up_rk4 + + # Make a combined u with the SBP and the DG parts. + u_comb = np.zeros(int(n_sbp_x*n_sbp_y) + nnodes) + u_comb[0:int(n_sbp_x*n_sbp_y)] = u_sbp + u_flat = flatten(u) + for i in range(int(n_sbp_x*n_sbp_y), int(n_sbp_x*n_sbp_y) + nnodes): + u_comb[i] = u_flat[i - int(n_sbp_x*n_sbp_y)].get() + dt_stepper = set_up_rk4("u", dt, u_comb, rhs) + + from grudge.shortcuts import make_visualizer + vis = make_visualizer(dcoll, vis_order=order) + + step = 0 + + # Create mesh for structured grid output + sbp_mesh = np.zeros((2, n_sbp_y, n_sbp_x)) + for j in range(0, n_sbp_y): + sbp_mesh[0, j, :] = x_sbp + for i in range(0, n_sbp_x): + sbp_mesh[1, :, i] = y_sbp + + for event in dt_stepper.run(t_end=final_time): + if isinstance(event, dt_stepper.StateComputed): + + step += 1 + + last_t = event.t + u_sbp = event.state_component[0:int(n_sbp_x*n_sbp_y)] + u_dg = event.state_component[int(n_sbp_x*n_sbp_y):] + + error_l2_dg = op.norm( + dcoll, + u_dg - u_analytic(nodes, t=last_t), + 2 + ) + print('DG L2 Error after Step ', step, error_l2_dg) + sbp_error = np.zeros((n_sbp_x*n_sbp_y)) + error_l2_sbp = 0 + for j in range(0, n_sbp_y): + for i in range(0, n_sbp_x): + sbp_error[i + j*n_sbp_x] = u_sbp[i + j*n_sbp_x].get() - \ + np.sin(10*(-c.dot([x_sbp[i], y_sbp[j]])/norm_c + + - last_t*norm_c)) + error_l2_sbp = error_l2_sbp + \ + dx*dy*(sbp_error[i + j*n_sbp_x]) ** 2 + + error_l2_sbp = np.sqrt(error_l2_sbp) + print('SBP L2 Error after Step ', step, error_l2_sbp) + + # Write out the DG data only + vis.write_vtk_file("dg-%04d.vtu" % step, + [("u", u_dg)], overwrite=True) + + # Try writing out a VTK file with the SBP data. + from pyvisfile.vtk import write_structured_grid + + # Overwrite existing files - this is annoying when debugging. + filename = "sbp_%04d.vts" % step + import os + if os.path.exists(filename): + os.remove(filename) + + write_structured_grid(filename, sbp_mesh, + point_data=[("u", u_sbp.get())]) + + +if __name__ == "__main__": + import sys + if len(sys.argv) > 1: + exec(sys.argv[1]) + else: + pytest.main([__file__]) From d886fcf3e0c6dc2589406704a45e349413d12119 Mon Sep 17 00:00:00 2001 From: Cory Mikida Date: Wed, 1 Sep 2021 15:18:35 -0500 Subject: [PATCH 2/7] Tests now run (but still fail) --- test/test_grudge.py | 35 ++++++++++++++++++++++------------- test/test_sbp_dg.py | 34 ++++++++++++++++++---------------- 2 files changed, 40 insertions(+), 29 deletions(-) diff --git a/test/test_grudge.py b/test/test_grudge.py index e72e30fa2..79e7555c5 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -980,7 +980,8 @@ def u_analytic(x, t=0): dudy_mat = -np.kron(d_y, np.eye(n_sbp_x)) # Number of nodes in our SBP-DG boundary discretization - sbp_nodes_y = thaw(sbp_bdry_discr.nodes(), actx)[1] + from meshmode.dof_array import flatten, unflatten + sbp_nodes_y = flatten(thaw(sbp_bdry_discr.nodes(), actx)[1]) # When projecting, we use nodes sorted in y, but we will have to unsort # afterwards to make sure projected solution is injected into DG BC # in the correct way. @@ -995,6 +996,7 @@ def u_analytic(x, t=0): dg_side_gg = np.zeros(int(west_nodes.shape[0]/(order+1))+1) counter = 0 for i in range(0, west_nodes.shape[0]): + west_nodes[i] = west_nodes[i].get() if i % (order+1) == 0: dg_side_gg[counter] = west_nodes[i] counter += 1 @@ -1041,7 +1043,7 @@ def rhs(t, u): dr_y = np.zeros(n_sbp_x*n_sbp_y) # Pull DG solution at western face to project. - u_dg_ts = u[int(n_sbp_x*n_sbp_y):].get() + u_dg_ts = u[int(n_sbp_x*n_sbp_y):] dg_west = np.zeros(nsbp_nodes) for i in range(0, nsbp_nodes): @@ -1083,11 +1085,7 @@ def rhs(t, u): # Add these at each point on the SBP half to get the SBP RHS. rhs_sbp = c[0]*dudx + c[1]*dudy - dl_x - dr_x - dl_y - dr_y - # Now pop this back into the device RHS vector. - rhs_sbp_dev = actx.np.zeros((n_sbp_x*n_sbp_y,)) - rhs_sbp_dev = rhs_sbp - - rhs_out[0:int(n_sbp_x*n_sbp_y)] = rhs_sbp_dev + rhs_out[0:int(n_sbp_x*n_sbp_y)] = rhs_sbp sbp_east = np.zeros(n_sbp_y) # Pull SBP domain values off of east face. @@ -1105,15 +1103,23 @@ def rhs(t, u): sbp_proj = sbp2dg.dot(sbp_east) # Second: Fix the ordering. sbp_proj = sbp_proj[unsort_args] + sbp_tag = dof_desc.DTAG_BOUNDARY("btag_sbp") + + u_dg_in = unflatten(actx, dcoll.discr_from_dd("vol"), + actx.from_numpy(u[int(n_sbp_x*n_sbp_y):])) + u_sbp_in = unflatten(actx, dcoll.discr_from_dd(sbp_tag), + actx.from_numpy(sbp_proj)) # Grudge DG RHS. # Critical step - now need to apply projected SBP state to the # proper nodal locations in u_dg. - rhs_out[int(n_sbp_x*n_sbp_y):] = adv_operator.operator( + dg_rhs = adv_operator.operator( t=t, - u=u[int(n_sbp_x*n_sbp_y):], - sbp_tag=sym.DTAG_BOUNDARY("btag_sbp"), std_tag=std_tag, - state_from_sbp=sbp_proj) + u=u_dg_in, + state_from_sbp=u_sbp_in, + sbp_tag=sbp_tag, std_tag=std_tag) + dg_rhs = flatten(dg_rhs) + rhs_out[int(n_sbp_x*n_sbp_y):] = dg_rhs.get() return rhs_out @@ -1123,8 +1129,9 @@ def rhs(t, u): # Make a combined u with the SBP and the DG parts. u_comb = actx.np.zeros(int(n_sbp_x*n_sbp_y) + nnodes) u_comb[0:int(n_sbp_x*n_sbp_y)] = u_sbp + u_flat = flatten(u) for i in range(int(n_sbp_x*n_sbp_y), int(n_sbp_x*n_sbp_y) + nnodes): - u_comb[i] = u[i - int(n_sbp_x*n_sbp_y)] + u_comb[i] = u_flat[i - int(n_sbp_x*n_sbp_y)].get() dt_stepper = set_up_rk4("u", dt, u_comb, rhs) from grudge.shortcuts import make_visualizer @@ -1169,9 +1176,11 @@ def rhs(t, u): # Write out the DG data if visualize: + u_dg_plot = unflatten(actx, dcoll.discr_from_dd("vol"), + actx.from_numpy(u_dg)) vis.write_vtk_file("eoc_dg-%s-%04d.vtu" % (nelem, step), - [("u", u_dg)], overwrite=True) + [("u", u_dg_plot)], overwrite=True) # Write out the SBP data. from pyvisfile.vtk import write_structured_grid diff --git a/test/test_sbp_dg.py b/test/test_sbp_dg.py index 09eec6245..d53e62dca 100644 --- a/test/test_sbp_dg.py +++ b/test/test_sbp_dg.py @@ -50,8 +50,6 @@ def test_sbp_dg(actx_factory, write_output=True, order=4): # DG Half. dim = 2 - from grudge import sym - nelem_x = 20 nelem_y = 20 mesh = mgen.generate_regular_rect_mesh(a=(0, -1), b=(1, 1), @@ -100,7 +98,6 @@ def u_analytic(x, t=0): from grudge.models.advection import WeakAdvectionSBPOperator std_tag = dof_desc.DTAG_BOUNDARY("btag_std") - from meshmode.mesh import BTAG_ALL adv_operator = WeakAdvectionSBPOperator(dcoll, c, inflow_u=lambda t: u_analytic( @@ -203,7 +200,7 @@ def u_analytic(x, t=0): dudy_mat = -np.kron(d_y, np.eye(n_sbp_x)) # Number of nodes in our SBP-DG boundary discretization - from meshmode.dof_array import flatten + from meshmode.dof_array import flatten, unflatten sbp_nodes_y = flatten(thaw(sbp_bdry_discr.nodes(), actx)[1]) # When projecting, we use nodes sorted in y, but we will have to unsort # afterwards to make sure projected solution is injected into DG BC @@ -283,11 +280,7 @@ def rhs(t, u): # Add these at each point on the SBP half to get the SBP RHS. rhs_sbp = c[0]*dudx + c[1]*dudy - dl_x - dr_x - dl_y - dr_y - # Now pop this back into the device RHS vector. - # rhs_sbp_dev = actx.np.zeros((n_sbp_x*n_sbp_y,)) - rhs_sbp_dev = rhs_sbp - - rhs_out[0:int(n_sbp_x*n_sbp_y)] = rhs_sbp_dev + rhs_out[0:int(n_sbp_x*n_sbp_y)] = rhs_sbp sbp_east = np.zeros(n_sbp_y) # Pull SBP domain values off of east face. @@ -305,15 +298,22 @@ def rhs(t, u): sbp_proj = sbp2dg.dot(sbp_east) # Second: Fix the ordering. sbp_proj = sbp_proj[unsort_args] + sbp_tag = dof_desc.DTAG_BOUNDARY("btag_sbp") + u_dg_in = unflatten(actx, dcoll.discr_from_dd("vol"), + actx.from_numpy(u[int(n_sbp_x*n_sbp_y):])) + u_sbp_in = unflatten(actx, dcoll.discr_from_dd(sbp_tag), + actx.from_numpy(sbp_proj)) # Grudge DG RHS. # Critical step - now need to apply projected SBP state to the # proper nodal locations in u_dg. - rhs_out[int(n_sbp_x*n_sbp_y):] = adv_operator.operator( + dg_rhs = adv_operator.operator( t=t, - u=u[int(n_sbp_x*n_sbp_y):], - state_from_sbp=sbp_proj, - sbp_tag=dof_desc.DTAG_BOUNDARY("btag_sbp"), std_tag=std_tag) + u=u_dg_in, + state_from_sbp=u_sbp_in, + sbp_tag=sbp_tag, std_tag=std_tag) + dg_rhs = flatten(dg_rhs) + rhs_out[int(n_sbp_x*n_sbp_y):] = dg_rhs.get() return rhs_out @@ -359,7 +359,7 @@ def rhs(t, u): error_l2_sbp = 0 for j in range(0, n_sbp_y): for i in range(0, n_sbp_x): - sbp_error[i + j*n_sbp_x] = u_sbp[i + j*n_sbp_x].get() - \ + sbp_error[i + j*n_sbp_x] = u_sbp[i + j*n_sbp_x] - \ np.sin(10*(-c.dot([x_sbp[i], y_sbp[j]])/norm_c + - last_t*norm_c)) error_l2_sbp = error_l2_sbp + \ @@ -369,8 +369,10 @@ def rhs(t, u): print('SBP L2 Error after Step ', step, error_l2_sbp) # Write out the DG data only + u_dg_plot = unflatten(actx, dcoll.discr_from_dd("vol"), + actx.from_numpy(u_dg)) vis.write_vtk_file("dg-%04d.vtu" % step, - [("u", u_dg)], overwrite=True) + [("u", u_dg_plot)], overwrite=True) # Try writing out a VTK file with the SBP data. from pyvisfile.vtk import write_structured_grid @@ -382,7 +384,7 @@ def rhs(t, u): os.remove(filename) write_structured_grid(filename, sbp_mesh, - point_data=[("u", u_sbp.get())]) + point_data=[("u", u_sbp)]) if __name__ == "__main__": From 2a241ca9004f19d739a450969b554371999f5b7a Mon Sep 17 00:00:00 2001 From: Cory Mikida Date: Wed, 1 Sep 2021 15:23:56 -0500 Subject: [PATCH 3/7] Minor changes to get the test to actually run --- test/test_grudge.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/test/test_grudge.py b/test/test_grudge.py index 79e7555c5..4f90fae96 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -838,7 +838,6 @@ def test_convergence_sbp_advec(actx_factory, order, order_sbp, spacing_factor, c from pytools.convergence import EOCRecorder eoc_rec = EOCRecorder() - from grudge import sym from sbp_operators import (sbp21, sbp42, sbp63) from projection import sbp_dg_projection @@ -917,7 +916,7 @@ def u_analytic(x, t=0): # Set up solution vector: # For now, timestep is the same as DG. - u_sbp = actx.np.zeros(int(n_sbp_x*n_sbp_y)) + u_sbp = np.zeros(int(n_sbp_x*n_sbp_y)) # Initial condition for j in range(0, n_sbp_y): @@ -1023,7 +1022,7 @@ def u_analytic(x, t=0): def rhs(t, u): # Initialize the entire RHS to 0. - rhs_out = actx.np.zeros(int(n_sbp_x*n_sbp_y) + int(nnodes)) + rhs_out = np.zeros(int(n_sbp_x*n_sbp_y) + int(nnodes)) # Fill the first part with the SBP half of the domain. @@ -1127,7 +1126,7 @@ def rhs(t, u): from grudge.shortcuts import set_up_rk4 # Make a combined u with the SBP and the DG parts. - u_comb = actx.np.zeros(int(n_sbp_x*n_sbp_y) + nnodes) + u_comb = np.zeros(int(n_sbp_x*n_sbp_y) + nnodes) u_comb[0:int(n_sbp_x*n_sbp_y)] = u_sbp u_flat = flatten(u) for i in range(int(n_sbp_x*n_sbp_y), int(n_sbp_x*n_sbp_y) + nnodes): From b7f7ff9625d202a061802d75110156554d018ed5 Mon Sep 17 00:00:00 2001 From: Cory Mikida Date: Thu, 2 Sep 2021 14:40:42 -0500 Subject: [PATCH 4/7] Minor test fixes --- test/test_grudge.py | 18 +++++++++--------- test/test_sbp_dg.py | 2 +- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/test/test_grudge.py b/test/test_grudge.py index 4f90fae96..cb84963b4 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -984,12 +984,12 @@ def u_analytic(x, t=0): # When projecting, we use nodes sorted in y, but we will have to unsort # afterwards to make sure projected solution is injected into DG BC # in the correct way. - nodesort = np.argsort(sbp_nodes_y.get()) + nodesort = np.argsort(sbp_nodes_y) nodesortlist = nodesort.tolist() - rangex = np.array(range(sbp_nodes_y.get().shape[0])) + rangex = np.array(range(sbp_nodes_y.shape[0])) unsort_args = [nodesortlist.index(x) for x in rangex] - west_nodes = np.sort(np.array(sbp_nodes_y.get())) + west_nodes = np.sort(np.array(sbp_nodes_y)) # Make element-aligned glue grid. dg_side_gg = np.zeros(int(west_nodes.shape[0]/(order+1))+1) @@ -1007,14 +1007,14 @@ def u_analytic(x, t=0): west_nodes) # Get mapping for western face - base_nodes = thaw(dcoll._volume_discr.nodes()) - nsbp_nodes = sbp_bdry_discr.nnodes + base_nodes = thaw(dcoll._volume_discr.nodes(), actx) + nsbp_nodes = sbp_bdry_discr.ndofs nodes_per_element = mesh.groups[0].nodes.shape[2] west_indices = np.zeros(nsbp_nodes) count = 0 # Sweep through first block of indices in the box. for i in range(0, (nelem-1)*2*nodes_per_element): - if base_nodes[0][i] < 1e-10: + if base_nodes[0][0][i][0] < 1e-10: # Make sure we're actually at the edge faces. if i % (2*nodes_per_element) < nodes_per_element: west_indices[count] = i @@ -1027,7 +1027,7 @@ def rhs(t, u): # Fill the first part with the SBP half of the domain. # Pull the SBP vector out of device array for now. - u_sbp_ts = u[0:int(n_sbp_x*n_sbp_y)].get() + u_sbp_ts = u[0:int(n_sbp_x*n_sbp_y)] dudx = np.zeros((n_sbp_x*n_sbp_y)) dudy = np.zeros((n_sbp_x*n_sbp_y)) @@ -1165,7 +1165,7 @@ def rhs(t, u): for j in range(0, n_sbp_y): for i in range(0, n_sbp_x): sbp_error[i + j*n_sbp_x] = \ - u_sbp[i + j*n_sbp_x].get() - \ + u_sbp[i + j*n_sbp_x] - \ np.sin(10*(-c.dot([x_sbp[i], y_sbp[j]]) / norm_c + last_t*norm_c)) error_l2_sbp = error_l2_sbp + \ @@ -1186,7 +1186,7 @@ def rhs(t, u): if visualize: filename = "eoc_sbp_%s-%04d.vts" % (nelem, step) write_structured_grid(filename, sbp_mesh, - point_data=[("u", u_sbp.get())]) + point_data=[("u", u_sbp)]) if c[0] > 0: eoc_rec.add_data_point(h, error_l2_dg) diff --git a/test/test_sbp_dg.py b/test/test_sbp_dg.py index d53e62dca..7d59af496 100644 --- a/test/test_sbp_dg.py +++ b/test/test_sbp_dg.py @@ -361,7 +361,7 @@ def rhs(t, u): for i in range(0, n_sbp_x): sbp_error[i + j*n_sbp_x] = u_sbp[i + j*n_sbp_x] - \ np.sin(10*(-c.dot([x_sbp[i], y_sbp[j]])/norm_c + - - last_t*norm_c)) + last_t*norm_c)) error_l2_sbp = error_l2_sbp + \ dx*dy*(sbp_error[i + j*n_sbp_x]) ** 2 From 1d5d342324648cb94956461638d696d99860adc4 Mon Sep 17 00:00:00 2001 From: Cory Mikida Date: Thu, 2 Sep 2021 14:49:31 -0500 Subject: [PATCH 5/7] Placate Flake8 --- test/test_sbp_dg.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_sbp_dg.py b/test/test_sbp_dg.py index 7d59af496..44c724e5b 100644 --- a/test/test_sbp_dg.py +++ b/test/test_sbp_dg.py @@ -360,8 +360,8 @@ def rhs(t, u): for j in range(0, n_sbp_y): for i in range(0, n_sbp_x): sbp_error[i + j*n_sbp_x] = u_sbp[i + j*n_sbp_x] - \ - np.sin(10*(-c.dot([x_sbp[i], y_sbp[j]])/norm_c + - last_t*norm_c)) + np.sin(10*(-c.dot([x_sbp[i], y_sbp[j]])/norm_c + + last_t*norm_c)) error_l2_sbp = error_l2_sbp + \ dx*dy*(sbp_error[i + j*n_sbp_x]) ** 2 From ae3d2193f93a2dea4b8fefe91a8fa58f3d7d754d Mon Sep 17 00:00:00 2001 From: Cory Mikida Date: Thu, 2 Sep 2021 15:51:00 -0500 Subject: [PATCH 6/7] Tests now pass! --- test/test_grudge.py | 14 +++++++------- test/test_sbp_dg.py | 7 +++---- 2 files changed, 10 insertions(+), 11 deletions(-) diff --git a/test/test_grudge.py b/test/test_grudge.py index cb84963b4..38ddf90b6 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -1007,14 +1007,14 @@ def u_analytic(x, t=0): west_nodes) # Get mapping for western face - base_nodes = thaw(dcoll._volume_discr.nodes(), actx) + base_nodes = flatten(thaw(dcoll._volume_discr.nodes(), actx)[0]) nsbp_nodes = sbp_bdry_discr.ndofs nodes_per_element = mesh.groups[0].nodes.shape[2] west_indices = np.zeros(nsbp_nodes) count = 0 # Sweep through first block of indices in the box. for i in range(0, (nelem-1)*2*nodes_per_element): - if base_nodes[0][0][i][0] < 1e-10: + if base_nodes[i] < 1e-10: # Make sure we're actually at the edge faces. if i % (2*nodes_per_element) < nodes_per_element: west_indices[count] = i @@ -1152,7 +1152,9 @@ def rhs(t, u): last_t = event.t u_sbp = event.state_component[0:int(n_sbp_x*n_sbp_y)] - u_dg = event.state_component[int(n_sbp_x*n_sbp_y):] + u_dg = unflatten(actx, dcoll.discr_from_dd("vol"), + actx.from_numpy( + event.state_component[int(n_sbp_x*n_sbp_y):])) error_l2_dg = op.norm( dcoll, @@ -1175,11 +1177,9 @@ def rhs(t, u): # Write out the DG data if visualize: - u_dg_plot = unflatten(actx, dcoll.discr_from_dd("vol"), - actx.from_numpy(u_dg)) vis.write_vtk_file("eoc_dg-%s-%04d.vtu" % (nelem, step), - [("u", u_dg_plot)], overwrite=True) + [("u", u_dg)], overwrite=True) # Write out the SBP data. from pyvisfile.vtk import write_structured_grid @@ -1189,7 +1189,7 @@ def rhs(t, u): point_data=[("u", u_sbp)]) if c[0] > 0: - eoc_rec.add_data_point(h, error_l2_dg) + eoc_rec.add_data_point(h, error_l2_dg.get()) else: eoc_rec.add_data_point(h, error_l2_sbp) diff --git a/test/test_sbp_dg.py b/test/test_sbp_dg.py index 44c724e5b..eb3f992e6 100644 --- a/test/test_sbp_dg.py +++ b/test/test_sbp_dg.py @@ -347,7 +347,8 @@ def rhs(t, u): last_t = event.t u_sbp = event.state_component[0:int(n_sbp_x*n_sbp_y)] - u_dg = event.state_component[int(n_sbp_x*n_sbp_y):] + u_dg = unflatten(actx, dcoll.discr_from_dd("vol"), + actx.from_numpy(event.state_component[int(n_sbp_x*n_sbp_y):])) error_l2_dg = op.norm( dcoll, @@ -369,10 +370,8 @@ def rhs(t, u): print('SBP L2 Error after Step ', step, error_l2_sbp) # Write out the DG data only - u_dg_plot = unflatten(actx, dcoll.discr_from_dd("vol"), - actx.from_numpy(u_dg)) vis.write_vtk_file("dg-%04d.vtu" % step, - [("u", u_dg_plot)], overwrite=True) + [("u", u_dg)], overwrite=True) # Try writing out a VTK file with the SBP data. from pyvisfile.vtk import write_structured_grid From f5abb2eea26ad99f4105d1f72dce9bb66ea795b5 Mon Sep 17 00:00:00 2001 From: Cory Mikida Date: Thu, 2 Sep 2021 18:02:45 -0500 Subject: [PATCH 7/7] Single quotes -> double quotes --- test/projection.py | 18 +++++++++--------- test/test_sbp_dg.py | 4 ++-- 2 files changed, 11 insertions(+), 11 deletions(-) diff --git a/test/projection.py b/test/projection.py index c3e1b5a59..c0bad137a 100644 --- a/test/projection.py +++ b/test/projection.py @@ -105,7 +105,7 @@ def make_projection(n, order): m = 2 # pb = 1 # from opt files - q = np.loadtxt('q_2.txt') + q = np.loadtxt("q_2.txt") r = 1 elif order == 4: # In K+W code, FD op size is n+1 @@ -113,7 +113,7 @@ def make_projection(n, order): m = 4 # pb = 2 # from opt files - q = np.loadtxt('q_4.txt') + q = np.loadtxt("q_4.txt") r = 5 elif order == 6: # In K+W code, FD op size is n+1 @@ -121,7 +121,7 @@ def make_projection(n, order): m = 6 # pb = 2 # from opt files - q = np.loadtxt('q_6.txt') + q = np.loadtxt("q_6.txt") r = 8 h = p_sbp @@ -147,7 +147,7 @@ def make_projection(n, order): pg2f[int(i[k])-1, int(j[k])-1] = qi[k] # Boundary solution. - qb = np.reshape(q[m*(p+1):], (r*(p+1), s,), order='F') + qb = np.reshape(q[m*(p+1):], (r*(p+1), s,), order="F") qb = np.transpose(qb) # Left block. @@ -155,8 +155,8 @@ def make_projection(n, order): qb = np.reshape((np.diag(2*np.mod(np.linspace(1, p+1, p+1, endpoint=True), 2) - 1)).dot(np.flipud(np.reshape(np.rot90(qb, - 2).transpose(), (p+1, r*s,), order='F'))), - (r*(p+1), s,), order='F').transpose() + 2).transpose(), (p+1, r*s,), order="F"))), + (r*(p+1), s,), order="F").transpose() # pg2f[np.arange(pg2f.shape[0]-s, pg2f.shape[0], 1), # np.arange(pg2f.shape[1]+1-r*(p+1)-1, pg2f.shape[1], 1)] = qb @@ -349,20 +349,20 @@ def sbp_sbp_test(): # we don't have those diagonal SBP operators coded yet. for qa in range(2, 6, 2): for qb in range(2, 6, 2): - print('Creating projection for (qa, qb) ', qa, qb) + print("Creating projection for (qa, qb) ", qa, qb) [pa2b, pb2a] = sbp_sbp_projection(na, qa, nb, qb) # np.savetxt('pa2b_test.txt', pa2b) # np.savetxt('pb2a_test.txt', pb2a) - print('Testing projection for (qa, qb) ', qa, qb) + print("Testing projection for (qa, qb) ", qa, qb) # Check the full operator for n in range(0, int(min(qa, qb)/2)-1): np.testing.assert_array_less(np.abs(pa2b.dot((xa ** n)) - xb ** n), 1000*eps) np.testing.assert_array_less(np.abs(pb2a.dot((xb ** n)) - xa ** n), 1000*eps) - print('Test passed for (qa, qb) ', qa, qb) + print("Test passed for (qa, qb) ", qa, qb) # Check the interior operator n_int = int(min(qa, qb))-1 diff --git a/test/test_sbp_dg.py b/test/test_sbp_dg.py index eb3f992e6..755a9ce33 100644 --- a/test/test_sbp_dg.py +++ b/test/test_sbp_dg.py @@ -355,7 +355,7 @@ def rhs(t, u): u_dg - u_analytic(nodes, t=last_t), 2 ) - print('DG L2 Error after Step ', step, error_l2_dg) + print("DG L2 Error after Step ", step, error_l2_dg) sbp_error = np.zeros((n_sbp_x*n_sbp_y)) error_l2_sbp = 0 for j in range(0, n_sbp_y): @@ -367,7 +367,7 @@ def rhs(t, u): dx*dy*(sbp_error[i + j*n_sbp_x]) ** 2 error_l2_sbp = np.sqrt(error_l2_sbp) - print('SBP L2 Error after Step ', step, error_l2_sbp) + print("SBP L2 Error after Step ", step, error_l2_sbp) # Write out the DG data only vis.write_vtk_file("dg-%04d.vtu" % step,