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

Solve DGL with non-constant (discrete) variable? #86

Open
stefmech opened this issue Jan 17, 2022 · 5 comments
Open

Solve DGL with non-constant (discrete) variable? #86

stefmech opened this issue Jan 17, 2022 · 5 comments

Comments

@stefmech
Copy link

Hi Mikael,

i would like to solve a DGL of this form:
.

is only defined by discrete values, so it cannot be incorparated as a expression. This is maybe a stupid question, but in the documentation I couldn't find a proper solution for such a problem. Is it possible to solve this kind of equation in shenfun with non-constant coefficients?

Best regards,
Steffen

@mikaem
Copy link
Member

mikaem commented Jan 18, 2022

Hi
You would be the first to implement something like this in shenfun, but it is not necessarily very difficult. You could probably use something like

and the inner product would then be

With test function and trial we get

So you get a dense matrix for each i. Perhaps not very clever after all?

@stefmech
Copy link
Author

Thank you for reply! Sorry for the not clearly formulated question. I mean, the non-constant variable is defined on the underlying grid N.

Would be the inner product not something like this?
.

Today, I saw that there's a scale functionality for an Expression, for example. What's about this? This could probably work for the above inner product, right?

Best regards

@mikaem
Copy link
Member

mikaem commented Jan 19, 2022

If q is defined on underlying grid, then there is nothing special to it. You can get it in spectral space simply by a forward transform. However, you do have something that needs to be computed like a nonlinear term. So I think you should compute it with something like

u = u_hat.backward() # get u on grid
qu = q*u # assuming u and q are now on the same grid in physical space
inner(v, qu)  

The scale functionality can only be used for a constant scaling of the entire expression. The scale can not be a function of x.

@stefmech
Copy link
Author

First of all, i would like to say thank you for your help! I saw your example "Demo - Lid driven cavity", where the nonlinear rhs was handled it in a similiar way as you mentioned above. However, I think the problem in your approach is , that the would result in a linear form, but not in a bilinear form. Or should I use an conjugate gradient solver which can use the matrix in the form "A*x"?
So would it be necessary to add a new case to the "inner"-function, e.g. setting up a bilinear form for the specific case, that there are two expressions consisting of a function and a testfunction?

I add an small example:

from shenfun import FunctionSpace, Function, grad, \
BlockMatrix, extract_bc_matrices, TestFunction, TrialFunction, inner, Array
import numpy as np

l = 100.
h = l/2
domain_x = (0, l)

N = 18
T = FunctionSpace(N, family='C', bc=(0,2), domain=domain_x)

u = TrialFunction(T)
v = TestFunction(T)

q = 100 # Is q is constant it can be simply multiplied

A = inner(q*grad(u), grad(v))

u_hat = Function(T)

# get boundary matrices
bc_mats = extract_bc_matrices([A])

# BlockMatrix for inhomogeneous part
BM = BlockMatrix(bc_mats)

# inhomogeneous part of solution
uh_hat = Function(T).set_boundary_dofs()

# additional part to be passed to the right hand side
b = Function(T)
# negative because added to right hand side
b = BM.matvec(-uh_hat, b)

# homogeneous part of solution
u_hat = A[0].solve(b)

u_ = u_hat.backward(kind='uniform')

print('Done')

@mikaem
Copy link
Member

mikaem commented Jan 20, 2022

Hi
You are right. Any nonlinear term or similar with multiplication of two functions need to use a linear form and go to the right hand side. So you would need to solve it iteratively. This is not just a shenfun problem, it is because when you multiply two spectral functions you need to compute a convolution and the order of the discretization goes naturally from N to 2N.

Dealiasing and pseudospectral techniques are usually the preferred way around this. Multiply u and q in real space (usually with 3/2N quadrature points) and then forward the outcome. There is no way I am aware of that can compute this easily in a bilinear form. In fact, it would take a trilinear form.

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