Forum Message

We have moved the forum to https://forum.ngsolve.org . This is an archived version of the topics until 05/05/23. All the topics were moved to the new forum and conversations can be continued there. This forum is just kept as legacy to not invalidate old links. If you want to continue a conversation just look for the topic in the new forum.

Notice

The forum is in read only mode.

Error L^2 norm for FESpace([V,F])

Amad
Topic Author
Offline
New Member

4 years 6 months ago - 4 years 6 months ago #1953 by Amad

Error L^2 norm for FESpace([V,F]) was created by Amad

Dear,

I've been working with NGSolve solving time harmonic elastic wave problem, with Robin boundary condition.
I am solving this problem using classical Galerkin and Hybrid Discontinuous Galerkin.

For the Galerkin method, it is working fine and the error in L2 norm is around 10^(-5).
For the HDG method, the solution seems correct (visually similar to the Galerkin method) but the L2 norm is not convergent.

For the HDG method, I choose the following spaces

Code:

V = VectorH1(mesh, order=order, complex=True)
F = VectorH1(mesh, order=order, complex=True)
fes = FESpace([V,F])

To compute the L2 norm, I am using the following code

Code:

u = CoefficientFunction(gfu.components[0])

#Error analysis
print ("Displacement L2-error:", sqrt (Integrate ( (u-uexact)*Conj(u-uexact), mesh)))

If I try to compute the norm using the following code, I get an error ("NgException: CompoundFESpace does not have an evaluator for VOL!")

Code:

gmesh = GridFunction(fes)
gmesh.Set(uexact)
u_exact = CoefficientFunction(gmesh)

u = CoefficientFunction(gfu.components[0])

print ("Displacement L2-error:", sqrt (Integrate ( (u-u_exact)*Conj(u-u_exact), mesh)))

How to compute the error in L^2 norm in fes = FESpace([V,F]) space?

Last edit: 4 years 6 months ago by Amad.

mneunteufel
Offline
Premium Member

4 years 6 months ago #1954 by mneunteufel

Replied by mneunteufel on topic Error L^2 norm for FESpace([V,F])

Dear Amad,

it is not possible to Set a Compound FESpace with one line

Code:

gmesh = GridFunction(fes)
gmesh.Set(uexact)

Setting each component individually

Code:

gmesh = GridFunction(fes)
gmesh.components[0].Set(uexact)

u = gfu.components[0]

print ("Displacement L2-error:", sqrt (Integrate ( (u-u_exact)*Conj(u-u_exact), mesh)))

should work.

Are you sure that for HDG your Facetvariable should also live in a VectorH1 space?

Code:

V = VectorH1(mesh, order=order, complex=True)
F = VectorH1(mesh, order=order, complex=True)
fes = FESpace([V,F])

Normally you want to have the dofs only on the facets like described here:
HDG tutorial

Maybe you'll need to make a Compound FESpace of three FacetFeSpaces.

Best,
Michael

The following user(s) said Thank You: Amad

Amad
Topic Author
Offline
New Member

4 years 6 months ago - 4 years 6 months ago #1958 by Amad

Replied by Amad on topic Error L^2 norm for FESpace([V,F])

Hi Michael,

Thank you for the prompt reply.
You are right, the chosen space was wrong.
But, I am facing a problem with the formulation.
As I am solving a vector problem, I need to choose a vector space. If I choose

Code:

V = VectorL2(mesh, order=order, complex=True)
F = VectorL2(mesh, order=order, complex=True)
fes = FESpace([V,F],dgjumps=True)

For the bilinear form

Code:

def epsilon(u):
    return 0.5 * (u.Deriv().trans + u.Deriv())
def sigma(u):
    return lamb*Trace(epsilon(u))*Id(mesh.dim) + 2*mu*epsilon(u)

a = BilinearForm(fes, condense=condense)
dS = dx(element_boundary=True)
a += InnerProduct(sigma(u), epsilon(v)) * dx - rho * omega**2 * InnerProduct(u, v) * dx +\
     beta*InnerProduct(jump_u,jump_v) * dS  +\
    -InnerProduct(sigma(u)*n, jump_v) * dS - InnerProduct(sigma(v)*n, jump_u) * dS
a += -1j*matA*uhat.Trace()*vhat.Trace() * ds

Then, I got the following error: "NgException: don't have a trace, primal evaluator = FN5ngfem6DiffOpINS_16DiffOpIdVectorH1ILi2ELNS_4VorBE0EEEEEvE"

If I choose the following spaces

Code:

V = VectorL2(mesh, order=order, complex=True)
F = FacetFESpace(mesh, order=order, complex=True)
fes = FESpace([V,F],dgjumps=True)

jump_u = u-uhat
jump_v = v-vhat

I got the error: "NgException: Dimensions don't match, op = - dims1 = 0: 2, dims2 = 0: 1"

Is it possible to choose a vector space for FacetFESpace?

Thank you in advance.

Last edit: 4 years 6 months ago by Amad.

mneunteufel
Offline
Premium Member

4 years 6 months ago #1963 by mneunteufel

Replied by mneunteufel on topic Error L^2 norm for FESpace([V,F])

Dear Amad,

you need to define 2 FacetFESpaces in 2D as the FacetFESpace is just a scalar. This explains the exception.
You can "construct" a VectorFacetSpace in the following way:

Code:

V = VectorL2(mesh, order=order, complex=True)
F = FacetFESpace(mesh, order=order, complex=True)
fes = FESpace([V,F,F],dgjumps=True)

(u,uhatx,uhaty), (v,vhatx,vhaty) = fes.TnT()

uhat = CoeffiecientFunction( (uhatx,uhaty) )
vhat = CoeffiecientFunction( (vhatx,vhaty) )

jump_u = u-uhat
jump_v = v-vhat

Best,
Michael

The following user(s) said Thank You: Amad

Amad
Topic Author
Offline
New Member

4 years 6 months ago #1964 by Amad

Replied by Amad on topic Error L^2 norm for FESpace([V,F])

Dear Michael,

Thank you for the reply.
Now there is no dimension problem. But, in the boundary terms, I had a problem.
It not recognise the uhat.Trace(), vhat.Trace().

Here is the bilinear form

Code:

def epsilon(u):
    return 0.5 * (u.Deriv().trans + u.Deriv())
def sigma(u):
    return lamb*Trace(epsilon(u))*Id(mesh.dim) + 2*mu*epsilon(u)

a = BilinearForm(fes, condense=condense)
dS = dx(element_boundary=True)
a += InnerProduct(sigma(u), epsilon(v)) * dx - rho * omega**2 * InnerProduct(u, v) * dx +\
     beta*InnerProduct(jump_u,jump_v) * dS  +\
    -InnerProduct(sigma(u)*n, jump_v) * dS - InnerProduct(sigma(v)*n, jump_u) * dS
a += -1j*matA*uhat.Trace()*vhat.Trace() * ds

And below is the error I got:

Code:

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-22-195408db1c27> in <module>()
     26 dS = dx(element_boundary=True)
     27 a += InnerProduct(sigma(u), epsilon(v)) * dx - rho * omega**2 *InnerProduct(u, v) * dx +     beta*InnerProduct(jump_u,jump_v) * dS +     -InnerProduct(sigma(u)*n, jump_v) * dS - InnerProduct(sigma(v)*n, jump_u) * dS
---> 28 a += -1j*matA*uhat.Trace()*vhat.Trace() * ds
     29 
     30 f = LinearForm(fes)

AttributeError: 'ngsolve.fem.CoefficientFunction' object has no attribute 'Trace'

mneunteufel
Offline
Premium Member

4 years 6 months ago #1965 by mneunteufel

Replied by mneunteufel on topic Error L^2 norm for FESpace([V,F])

Hi Amad,

the exception explains that CoeffiientFunctions do not have a trace. Only test- and trialfunctions do have a trace operator defined. You need to construct it by yourself by

Code:

uhat = CoeffiecientFunction( (uhatx,uhaty) )
vhat = CoeffiecientFunction( (vhatx,vhaty) )

uhat_tr = CoeffiecientFunction( (uhatx.Trace(),uhaty.Trace()) )
vhat_tr = CoeffiecientFunction( (vhatx.Trace(),vhaty.Trace()) )

...

a += -1j*matA*uhat_tr*vhat_tr* ds

Best,
Michael

The following user(s) said Thank You: Amad

Time to create page: 0.183 seconds

Login

NGSolve on Twitter