What can I do to calculate the L2 norm between the exact solution and my numerical solution?

Specific Domains Numerics
, ,
1

I am solving a differential equation by the MOL method, it has an exact solution so I want to calculate the error in the L2 norm. I made this code

F1(x) = A*(cosh(B*x))
F2(y) = A*(cosh(B*y))
FC = -(2*(B^2)*K - W)/(ρ*C_s)    
F3(t::Real) = t > t_c ? E * exp(FC * t) : 1 + E * exp(t)
T(t::Real)= x -> F1(x[1])*F2(x[2])*F3(t)
(...)
uₕₜ = solve(ode_solver,op_Af,u₀,t_0,t_f)

  error_L2 = Vector{Float64}(undef, 22)
for (uₕ,t) in uₕₜ
      error_L2[i] = √(integrate((T(t) - uₕ)^2, dΩ))
end

but when trying to up 2 the difference between the exact solution and my approximate solution it gives me an error. What can I do to calculate this norm?

2

Hi,

welcome to the community :slight_smile:

It makes it easier for us to help, if you post a fully runnable snippet of your code, including dependencies.

And can you provide the error message too?

3

Thanks for answering, of course I’ll put the code:

using GeometryBasics
using   Gridap,
        Gridap.Algebra,
        Gridap.FESpaces,
        Gridap.ReferenceFEs,
        Gridap.Arrays,
        Gridap.Geometry,
        Gridap.Fields,
        Gridap.CellData,
        GridapGmsh
using IterativeSolvers
using IncompleteLU
using SparseArrays
using BenchmarkTools 
using Plots
using CPUTime
using LinearAlgebra

output_path = "output_file/"
t_f = 1.
t_c = 0.2
Δt = 0.05
ρ = 1000            
C_s = 4200        
K = 84030       

F1(x) = A*(cosh(B*x))
F2(y) = A*(cosh(B*y))
FC = -(2*(B^2)*K - W)/(ρ*C_s)    
F3(t::Real) = t > t_c ? E * exp(FC * t) : 1 + E * exp(t)  
Q(x, t::Real) = t > t_c ? 0 : -C_s * F1(x[1]) * F2(x[2])
T(t::Real)= x -> F1(x[1])*F2(x[2])*F3(t)
T_0(x, t::Real) = t > t_c ? (1 + E * exp(t_c)) * F1(x[1]) * F2(x[2]) : (1 + E) * F1(x[1]) * F2(x[2])
T_in(t::Real) = x -> T_0(x,t)
g(x, t::Real) = F1(x[1]) * F2(x[2]) * F3(t)
g(t::Real) = x -> g(x,t)

model = CartesianDiscreteModel((0,1,0,1),(20,20))
Ω = Interior(model)
order = 2;
reffe = ReferenceFE(lagrangian, Float64, order)
dΩ = Measure(Ω, 2*order)
V = TestFESpace(model, reffe, conformity = :H1, dirichlet_tags = ["Gamma_D", "Gamma_N", "Gamma_R"], vector_type = Vector{Float64})
U = TransientTrialFESpace(V,g)

Qsource(t::Real) = x -> Q(x, t)
K_bil(t,u,v) = ∫(K*(∇(u) ⋅ ∇(v)) + W*(u*v))*dΩ 
M_bil(t,u,v) = ∫(ρ*C_s*(u*v))*dΩ
b_lin(t,v) = ∫( Qsource(t)*v )*dΩ

op_Af = TransientAffineFEOperator(M_bil, K_bil, b_lin, U, V)
linear_solver = LUSolver()
ode_solver = ThetaMethod(linear_solver,Δt,1)

t_0 = 0.0
u₀ = interpolate(T_in(t_0),U(t_0))
uₕₜ = solve(ode_solver,op_Af,u₀,t_0,t_f)

error_L2 = Vector{Float64}(undef, 22)
for (uₕ,t) in uₕₜ
      error_L2[i] = √(integrate((T(t) - uₕ)^2, dΩ))
end

In the part of calculating the error in the L2 norm it gives me the following error:

ERROR: MethodError: no method matching ^(::Gridap.CellData.OperationCellField{ReferenceDomain}, ::Int64)
Closest candidates are:
^(::Union{AbstractChar, AbstractString}, ::Integer) at strings/basic.jl:718
^(::AbstractGray{T} where T, ::Integer) at C:\Users\Manuel.julia\packages\ColorVectorSpace\QI5vM\src\ColorVectorSpace.jl:333
^(::AbstractGray{T} where T, ::Real) at C:\Users\Manuel.julia\packages\ColorVectorSpace\QI5vM\src\ColorVectorSpace.jl:334

4

An example is given at the end of this tutorial.

I think you need something like sum(abs2(T - uₕ) * dΩ), though in the tutorial it first projects the analytical solution into the FEM basis. In general, if you want to compose an arbitrary function f with a cell field u, you need to use f ∘ u, e.g. (x -> x^2) ∘ uₕ. But a specialized CellField method is already defined for abs2.

1 Like