Mercurial > repos > public > sbplib_julia

comparison test/testSbpOperators.jl @ 678:730565f7cc2e feature/laplace_opset

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Update constructor tests for Laplace
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Sat, 06 Feb 2021 15:18:05 +0100
parents 538ccbaeb1f8
children 54ce3f9771e5
comparison
equal deleted inserted replaced
677:011863b3f24c 678:730565f7cc2e
334 g_1D = EquidistantGrid(101, 0.0, 1.) 334 g_1D = EquidistantGrid(101, 0.0, 1.)
335 g_3D = EquidistantGrid((51,101,52), (0.0, -1.0, 0.0), (1., 1., 1.)) 335 g_3D = EquidistantGrid((51,101,52), (0.0, -1.0, 0.0), (1., 1., 1.))
336 @testset "Constructors" begin 336 @testset "Constructors" begin
337 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) 337 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
338 @testset "1D" begin 338 @testset "1D" begin
339 L = Laplace(g_1D, op.innerStencil, op.closureStencils) 339 # Create all tensor mappings included in Laplace
340 Δ = laplace(g_1D, op.innerStencil, op.closureStencils)
341 H = quadrature(g_1D, op.quadratureClosure)
342 Hi = InverseDiagonalQuadrature(g_1D, op.quadratureClosure)
343
344 (id_l, id_r) = boundary_identifiers(g_1D)
345
346 e_l = BoundaryRestriction(g_1D,op.eClosure,id_l)
347 e_r = BoundaryRestriction(g_1D,op.eClosure,id_r)
348 e_dict = Dict(Pair(id_l,e_l),Pair(id_r,e_r))
349
350 d_l = NormalDerivative(g_1D,op.dClosure,id_l)
351 d_r = NormalDerivative(g_1D,op.dClosure,id_r)
352 d_dict = Dict(Pair(id_l,d_l),Pair(id_r,d_r))
353
354 H_l = boundary_quadrature(g_1D,op.quadratureClosure,id_r)
355 H_r = boundary_quadrature(g_1D,op.quadratureClosure,id_r)
356 Hb_dict = Dict(Pair(id_l,H_l),Pair(id_r,H_r))
357
358 # TODO: Not sure why this doesnt work? Comparing the fields of
359 # Laplace seems to work. Reformulate below once solved.
360 @test_broken Laplace(Δ,H,Hi,e_dict,d_dict,Hb_dict) == Laplace(g_1D, sbp_operators_path()*"standard_diagonal.toml"; order=4)
361 L = Laplace(g_1D, sbp_operators_path()*"standard_diagonal.toml"; order=4)
362 @test L.D == Δ
363 @test L.H == H
364 @test L.H_inv == Hi
365 @test L.e == e_dict
366 @test L.d == d_dict
367 @test L.H_boundary == Hb_dict
368
369 @test L isa TensorMapping{T,1,1} where T
370 end
371 @testset "3D" begin
372 # Create all tensor mappings included in Laplace
373 Δ = laplace(g_3D, op.innerStencil, op.closureStencils)
374 H = quadrature(g_3D, op.quadratureClosure)
375 Hi = InverseDiagonalQuadrature(g_3D, op.quadratureClosure)
376
377 (id_l, id_r, id_s, id_n, id_b, id_t) = boundary_identifiers(g_3D)
378
379 e_l = BoundaryRestriction(g_3D,op.eClosure,id_l)
380 e_r = BoundaryRestriction(g_3D,op.eClosure,id_r)
381 e_s = BoundaryRestriction(g_3D,op.eClosure,id_s)
382 e_n = BoundaryRestriction(g_3D,op.eClosure,id_n)
383 e_b = BoundaryRestriction(g_3D,op.eClosure,id_b)
384 e_t = BoundaryRestriction(g_3D,op.eClosure,id_t)
385 e_dict = Dict(Pair(id_l,e_l),Pair(id_r,e_r),
386 Pair(id_s,e_s),Pair(id_n,e_n),
387 Pair(id_b,e_b),Pair(id_t,e_t))
388
389 d_l = NormalDerivative(g_3D,op.dClosure,id_l)
390 d_r = NormalDerivative(g_3D,op.dClosure,id_r)
391 d_s = NormalDerivative(g_3D,op.dClosure,id_s)
392 d_n = NormalDerivative(g_3D,op.dClosure,id_n)
393 d_b = NormalDerivative(g_3D,op.dClosure,id_b)
394 d_t = NormalDerivative(g_3D,op.dClosure,id_t)
395 d_dict = Dict(Pair(id_l,d_l),Pair(id_r,d_r),
396 Pair(id_s,d_s),Pair(id_n,d_n),
397 Pair(id_b,d_b),Pair(id_t,d_t))
398
399 H_l = boundary_quadrature(g_3D,op.quadratureClosure,id_r)
400 H_r = boundary_quadrature(g_3D,op.quadratureClosure,id_r)
401 H_s = boundary_quadrature(g_3D,op.quadratureClosure,id_s)
402 H_n = boundary_quadrature(g_3D,op.quadratureClosure,id_n)
403 H_b = boundary_quadrature(g_3D,op.quadratureClosure,id_b)
404 H_t = boundary_quadrature(g_3D,op.quadratureClosure,id_t)
405 Hb_dict = Dict(Pair(id_l,H_l),Pair(id_r,H_r),
406 Pair(id_s,H_s),Pair(id_n,H_n),
407 Pair(id_b,H_b),Pair(id_t,H_t))
408
409 # TODO: Not sure why this doesnt work? Comparing the fields of
410 # Laplace seems to work. Reformulate below once solved.
411 @test_broken Laplace(Δ,H,Hi,e_dict,d_dict,Hb_dict) == Laplace(g_3D, sbp_operators_path()*"standard_diagonal.toml"; order=4)
412 L = Laplace(g_3D, sbp_operators_path()*"standard_diagonal.toml"; order=4)
413 @test L.D == Δ
414 @test L.H == H
415 @test L.H_inv == Hi
416 @test L.e == e_dict
417 @test L.d == d_dict
418 @test L.H_boundary == Hb_dict
419 @test L isa TensorMapping{T,3,3} where T
420 end
421 end
422
423 @testset "laplace" begin
424 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
425 @testset "1D" begin
426 L = laplace(g_1D, op.innerStencil, op.closureStencils)
340 @test L == SecondDerivative(g_1D, op.innerStencil, op.closureStencils) 427 @test L == SecondDerivative(g_1D, op.innerStencil, op.closureStencils)
341 @test L isa TensorMapping{T,1,1} where T 428 @test L isa TensorMapping{T,1,1} where T
342 end 429 end
343 @testset "3D" begin 430 @testset "3D" begin
344 L = Laplace(g_3D, op.innerStencil, op.closureStencils) 431 L = laplace(g_3D, op.innerStencil, op.closureStencils)
345 @test L isa TensorMapping{T,3,3} where T
346 Dxx = SecondDerivative(g_3D, op.innerStencil, op.closureStencils,1) 432 Dxx = SecondDerivative(g_3D, op.innerStencil, op.closureStencils,1)
347 Dyy = SecondDerivative(g_3D, op.innerStencil, op.closureStencils,2) 433 Dyy = SecondDerivative(g_3D, op.innerStencil, op.closureStencils,2)
348 Dzz = SecondDerivative(g_3D, op.innerStencil, op.closureStencils,3) 434 Dzz = SecondDerivative(g_3D, op.innerStencil, op.closureStencils,3)
349 @test L == Dxx + Dyy + Dzz 435 @test L == Dxx + Dyy + Dzz
350 end 436 @test L isa TensorMapping{T,3,3} where T
437 end
438 end
439
440 @testset "quadrature" begin
441 end
442
443 @testset "inverse_quadrature" begin
444 end
445
446 @testset "boundary_restriction" begin
447 end
448
449 @testset "normal_restriction" begin
450 end
451
452 @testset "boundary_quadrature" begin
351 end 453 end
352 454
353 # Exact differentiation is measured point-wise. In other cases 455 # Exact differentiation is measured point-wise. In other cases
354 # the error is measured in the l2-norm. 456 # the error is measured in the l2-norm.
355 @testset "Accuracy" begin 457 @testset "Accuracy" begin
364 Δv = evalOn(g_3D,(x,y,z) -> -sin(x) - cos(y) + exp(z)) 466 Δv = evalOn(g_3D,(x,y,z) -> -sin(x) - cos(y) + exp(z))
365 467
366 # 2nd order interior stencil, 1st order boundary stencil, 468 # 2nd order interior stencil, 1st order boundary stencil,
367 # implies that L*v should be exact for binomials up to order 2. 469 # implies that L*v should be exact for binomials up to order 2.
368 @testset "2nd order" begin 470 @testset "2nd order" begin
369 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2) 471 L = Laplace(g_3D,sbp_operators_path()*"standard_diagonal.toml"; order=2)
370 L = Laplace(g_3D,op.innerStencil,op.closureStencils)
371 @test L*polynomials[1] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9 472 @test L*polynomials[1] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
372 @test L*polynomials[2] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9 473 @test L*polynomials[2] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
373 @test L*polynomials[3] ≈ polynomials[1] atol = 5e-9 474 @test L*polynomials[3] ≈ polynomials[1] atol = 5e-9
374 @test L*v ≈ Δv rtol = 5e-2 norm = l2 475 @test L*v ≈ Δv rtol = 5e-2 norm = l2
375 end 476 end
376 477
377 # 4th order interior stencil, 2nd order boundary stencil, 478 # 4th order interior stencil, 2nd order boundary stencil,
378 # implies that L*v should be exact for binomials up to order 3. 479 # implies that L*v should be exact for binomials up to order 3.
379 @testset "4th order" begin 480 @testset "4th order" begin
380 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) 481 L = Laplace(g_3D,sbp_operators_path()*"standard_diagonal.toml"; order=4)
381 L = Laplace(g_3D,op.innerStencil,op.closureStencils)
382 # NOTE: high tolerances for checking the "exact" differentiation 482 # NOTE: high tolerances for checking the "exact" differentiation
383 # due to accumulation of round-off errors/cancellation errors? 483 # due to accumulation of round-off errors/cancellation errors?
384 @test L*polynomials[1] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9 484 @test L*polynomials[1] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
385 @test L*polynomials[2] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9 485 @test L*polynomials[2] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
386 @test L*polynomials[3] ≈ polynomials[1] atol = 5e-9 486 @test L*polynomials[3] ≈ polynomials[1] atol = 5e-9