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

First Order Schemes #64

Open
jettrich opened this issue Sep 4, 2022 · 1 comment
Open

First Order Schemes #64

jettrich opened this issue Sep 4, 2022 · 1 comment

Comments

@jettrich
Copy link

jettrich commented Sep 4, 2022

Hi,
in case I want to use findiff for educational purpose, showing the differences in accuracy for first, second and higher order schemes, how can I enforce the usage of pure forward/backward differencing? Furthermore, e.g. thinking of advection, is it possible to realize upwind or TVD schemes with findiff?
Any help is greatly appreciated, kind regards!

@maroba
Copy link
Owner

maroba commented Sep 5, 2022

Hi,

if you want to enforce a specific difference scheme, you can use the Stencil class in findiff.stencils. Assuming you are using the currently latest version, 0.9.2, an example in 1D would read:

import numpy as np
from findiff.stencils import Stencil

# some example data:
x = np.linspace(0, 1, 101)
dx = x[1] - x[0]
f = np.sin(x)

# define an explicit forward stencil for the first derivative along axis 0:
stencil = Stencil(offsets=[0, 1, 2], partials:{(1,): 1}, spacings=dx)

Note that the partials keyword is kind of awkward. It allows you to define
a linear combination of partial derivative with constant coefficients with one single dict. The keys are tuples and each item of the tuple specifies the degree of the derivative along the axis. The values are the constant coefficients. That is, if you have a 2D grid, a partial derivative derivative along the axis 1 would read partials={(0, 1): 1} and a mixed one (1, 1): 1. Anyways, continuing with the 1D example:

stencil.accuracy
Out: 2

and

stencil.values
Out: {(0,): -150.0, (1,): 200.0, (2,): -50.0}

Apply the stencil at a given point, e.g. at x[5]:

stencil(f, at=5)
Out: 0.9987835384105308

Doing this point-wise is inefficient, though. But you can also apply the stencil to slices and masks instead of single points by using the on argument instead of at. This is much faster, of course.

I guess one could easily realize an upwind or TVD schemes using the Stencil class.

By the way, I'm about to release a new major version (1.0.0) in the next two weeks or so. It will contain symbolic capabilities which allows one to derive discretization schemes and analyze their stability behavior analytically. Maybe that's interesting for you.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants