kinnala / scikit-fem

Simple finite element assemblers
https://scikit-fem.readthedocs.io
BSD 3-Clause "New" or "Revised" License
488 stars 79 forks source link

advection problem #649

Closed nschloe closed 3 years ago

nschloe commented 3 years ago

I tried to formulate a diffusion-advection problem

-\Delta u + a\cdot \nabla u = f

(math rendered with Purple Pi) in skfem:

import meshzoo
from skfem import BilinearForm, LinearForm, MeshTri, InteriorBasis, ElementTriP1
from skfem.helpers import dot, grad

@BilinearForm
def laplace(u, v, _):
    return dot(grad(u), grad(v)) + dot([1.0, 1.0], u) * v

@LinearForm
def rhs(v, _):
    return 1.0 * v

points, cells = meshzoo.disk(6, 100)
m = MeshTri(points.T.copy(), cells.T.copy())
basis = InteriorBasis(m, ElementTriP1())
A = laplace.assemble(basis)
b = rhs.assemble(basis)
A, b = enforce(A, b, D=m.boundary_nodes())

This results in the cryptic error

ValueError: operands could not be broadcast together with remapped shapes [original->remapped]: (2,)->(2) (60000,3)->(3,60000)

Any hint on what's going wrong?

kinnala commented 3 years ago

I think dot([1., 1.], u) causes it. I'm quite sure dot-helper has no special case for such constant vectors. I'll try to consider how it could work.

kinnala commented 3 years ago

One question though. The test function u is scalar, right? So what do you mean by a dot u?

kinnala commented 3 years ago

Above I meant "trial function u".

Trying to recover from wisdom teeth removal :0).

nschloe commented 3 years ago

One question though. The test function u is scalar, right? So what do you mean by a dot u?

It should be

dot([1,1], grad(u))

and this works indeed. This one's on me, thanks for the quick reply!

gdmcbain commented 3 years ago

See also ex 25 Forced convection. This is just pure Galerkin. Longer term, I'm looking to investigate more specialized discretizations : method of characteristics #600, streamwise upwind Petrov Galerkin, discontinuous Galerkin, ...