Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Make operator for evaluating as ncc #253

Open
kburns opened this issue Jun 20, 2023 · 0 comments
Open

Make operator for evaluating as ncc #253

kburns opened this issue Jun 20, 2023 · 0 comments
Labels
enhancement New feature or request

Comments

@kburns
Copy link
Member

kburns commented Jun 20, 2023

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'])
@kburns kburns added the enhancement New feature or request label Jun 20, 2023
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement New feature or request
Projects
None yet
Development

No branches or pull requests

1 participant