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

Amad
Topic Author
Offline
New Member

5 years 1 month ago - 5 years 1 month 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: 5 years 1 month ago by Amad.

mneunteufel
Offline
Premium Member

5 years 1 month 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

5 years 1 month ago - 5 years 1 month 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: 5 years 1 month ago by Amad.

mneunteufel
Offline
Premium Member

5 years 1 month 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

5 years 1 month 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

5 years 1 month 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.105 seconds