It might be nice to create an operator / easy interface for evaluating multiplications using NCC matrices. We have this capability with the evaluate_as_ncc method, but not an operator that uses that. For instance, if you have a PDE LHS=RHS that satisfies some conservation law, this won't always be true if you multiple the equation through by an NCC as NCC*LHS=NCC*RHS. This is due to aliasing and truncation issues on the (a0,b0) Jacobi grid. Here's a test that gets the right result, but isn't pretty:
def test_solve_jacobi_ncc(N, a0, b0, k_ncc, k_arg, dealias, dtype):
"""Solve f(x) * u(x) = f(x) * g(x)."""
c, d, b = build_jacobi(N, a0, b0, (0, 1), dealias, dtype)
b_ncc = b.clone_with(a=a0+k_ncc, b=b0+k_ncc)
b_arg = b.clone_with(a=a0+k_arg, b=b0+k_arg)
f = d.Field(bases=b_ncc)
g = d.Field(bases=b_arg)
u = d.Field(bases=b_arg)
f.fill_random('g')
g.fill_random('g')
vars = [g]
w0 = f * g
w1 = w0.reinitialize(ncc=True, ncc_vars=vars)
problem = d3.LBVP(vars)
problem.add_equation((w1, 0))
solver = problem.build_solver()
w1.store_ncc_matrices(vars, solver.subproblems)
w0 = w0.evaluate()
w1 = w1.evaluate_as_ncc()
problem = d3.LBVP([u])
problem.add_equation((f*u, w1))
solver = problem.build_solver()
solver.solve()
assert np.allclose(u['c'], g['c'])
It might be nice to create an operator / easy interface for evaluating multiplications using NCC matrices. We have this capability with the
evaluate_as_nccmethod, but not an operator that uses that. For instance, if you have a PDELHS=RHSthat satisfies some conservation law, this won't always be true if you multiple the equation through by an NCC asNCC*LHS=NCC*RHS. This is due to aliasing and truncation issues on the (a0,b0) Jacobi grid. Here's a test that gets the right result, but isn't pretty: