|
| 1 | +module MaxwellNonConformingOctreeModelsTests |
| 2 | + using P4est_wrapper |
| 3 | + using GridapP4est |
| 4 | + using Gridap |
| 5 | + using PartitionedArrays |
| 6 | + using GridapDistributed |
| 7 | + using MPI |
| 8 | + using Gridap.FESpaces |
| 9 | + using FillArrays |
| 10 | + using Logging |
| 11 | + using Test |
| 12 | + |
| 13 | + function test_transfer_ops_and_redistribute(ranks, |
| 14 | + dmodel::GridapDistributed.DistributedDiscreteModel{Dc}, |
| 15 | + order) where Dc |
| 16 | + ref_coarse_flags=map(ranks,partition(get_cell_gids(dmodel.dmodel))) do rank,indices |
| 17 | + flags=zeros(Cint,length(indices)) |
| 18 | + flags.=nothing_flag |
| 19 | + |
| 20 | + flags[1]=refine_flag |
| 21 | + flags[own_length(indices)]=refine_flag |
| 22 | + |
| 23 | + # To create some unbalance |
| 24 | + if (rank%2==0 && own_length(indices)>1) |
| 25 | + flags[div(own_length(indices),2)]=refine_flag |
| 26 | + end |
| 27 | + flags |
| 28 | + end |
| 29 | + fmodel,glue=Gridap.Adaptivity.adapt(dmodel,ref_coarse_flags); |
| 30 | + |
| 31 | + # Solve coarse |
| 32 | + uH,UH=solve_maxwell(dmodel,order) |
| 33 | + check_error_maxwell(dmodel,order,uH) |
| 34 | + |
| 35 | + # Solve fine |
| 36 | + uh,Uh=solve_maxwell(fmodel,order) |
| 37 | + check_error_maxwell(fmodel,order,uh) |
| 38 | + |
| 39 | + Ωh = Triangulation(fmodel) |
| 40 | + degree = 2*(order+1) |
| 41 | + dΩh = Measure(Ωh,degree) |
| 42 | + |
| 43 | + # prolongation via interpolation |
| 44 | + uHh=interpolate(uH,Uh) |
| 45 | + e = uh - uHh |
| 46 | + el2 = sqrt(sum( ∫( e⋅e )*dΩh )) |
| 47 | + tol=1e-6 |
| 48 | + println("[INTERPOLATION] el2 < tol: $(el2) < $(tol)") |
| 49 | + @assert el2 < tol |
| 50 | + |
| 51 | + # prolongation via L2-projection |
| 52 | + # Coarse FEFunction -> Fine FEFunction, by projection |
| 53 | + ahp(u,v) = ∫(v⋅u)*dΩh |
| 54 | + lhp(v) = ∫(v⋅uH)*dΩh |
| 55 | + oph = AffineFEOperator(ahp,lhp,Uh,Uh) |
| 56 | + uHh = solve(oph) |
| 57 | + e = uh - uHh |
| 58 | + el2 = sqrt(sum( ∫( e⋅e )*dΩh )) |
| 59 | + println("[L2 PROJECTION] el2 < tol: $(el2) < $(tol)") |
| 60 | + @assert el2 < tol |
| 61 | + |
| 62 | + # restriction via interpolation |
| 63 | + uhH=interpolate(uh,UH) |
| 64 | + e = uH - uhH |
| 65 | + el2 = sqrt(sum( ∫( e⋅e )*dΩh )) |
| 66 | + println("[INTERPOLATION] el2 < tol: $(el2) < $(tol)") |
| 67 | + @assert el2 < tol |
| 68 | + |
| 69 | + # restriction via L2-projection |
| 70 | + ΩH = Triangulation(dmodel) |
| 71 | + degree = 2*(order+1) |
| 72 | + dΩH = Measure(ΩH,degree) |
| 73 | + |
| 74 | + dΩhH = Measure(ΩH,Ωh,2*order) |
| 75 | + aHp(u,v) = ∫(v⋅u)*dΩH |
| 76 | + lHp(v) = ∫(v⋅uh)*dΩhH |
| 77 | + oph = AffineFEOperator(aHp,lHp,UH,UH) |
| 78 | + uhH = solve(oph) |
| 79 | + e = uH - uhH |
| 80 | + el2 = sqrt(sum( ∫( e⋅e )*dΩH )) |
| 81 | + |
| 82 | + fmodel_red, red_glue=GridapDistributed.redistribute(fmodel); |
| 83 | + uh_red,Uh_red=solve_maxwell(fmodel_red,order) |
| 84 | + check_error_maxwell(fmodel_red,order,uh_red) |
| 85 | + |
| 86 | + trian = Triangulation(fmodel_red) |
| 87 | + degree = 2*(order+1) |
| 88 | + dΩhred = Measure(trian,degree) |
| 89 | + |
| 90 | + u_ex, f_ex=get_analytical_functions(Dc) |
| 91 | + |
| 92 | + uhred2 = GridapDistributed.redistribute_fe_function(uh,Uh_red,fmodel_red,red_glue) |
| 93 | + |
| 94 | + e = u_ex - uhred2 |
| 95 | + el2 = sqrt(sum( ∫( e⋅e )*dΩhred )) |
| 96 | + println("[REDISTRIBUTE SOLUTION] el2 < tol: $(el2) < $(tol)") |
| 97 | + @assert el2 < tol |
| 98 | + |
| 99 | + fmodel_red |
| 100 | + end |
| 101 | + |
| 102 | + function test_refine_and_coarsen_at_once(ranks, |
| 103 | + dmodel::OctreeDistributedDiscreteModel{Dc}, |
| 104 | + order) where Dc |
| 105 | + degree = 2*order+1 |
| 106 | + ref_coarse_flags=map(ranks,partition(get_cell_gids(dmodel.dmodel))) do rank,indices |
| 107 | + flags=zeros(Cint,length(indices)) |
| 108 | + flags.=nothing_flag |
| 109 | + if (rank==1) |
| 110 | + flags[1:min(2^Dc,own_length(indices))].=coarsen_flag |
| 111 | + end |
| 112 | + flags[own_length(indices)]=refine_flag |
| 113 | + flags |
| 114 | + end |
| 115 | + fmodel,glue=Gridap.Adaptivity.adapt(dmodel,ref_coarse_flags); |
| 116 | + |
| 117 | + # Solve coarse |
| 118 | + uH,UH=solve_maxwell(dmodel,order) |
| 119 | + check_error_maxwell(dmodel,order,uH) |
| 120 | + |
| 121 | + # Solve fine |
| 122 | + uh,Uh=solve_maxwell(fmodel,order) |
| 123 | + check_error_maxwell(fmodel,order,uh) |
| 124 | + |
| 125 | + # # prolongation via interpolation |
| 126 | + uHh=interpolate(uH,Uh) |
| 127 | + e = uh - uHh |
| 128 | + |
| 129 | + trian = Triangulation(fmodel) |
| 130 | + degree = 2*(order+1) |
| 131 | + dΩh = Measure(trian,degree) |
| 132 | + |
| 133 | + el2 = sqrt(sum( ∫( e⋅e )*dΩh )) |
| 134 | + tol=1e-6 |
| 135 | + @assert el2 < tol |
| 136 | + end |
| 137 | + |
| 138 | + function test_2d(ranks,order) |
| 139 | + coarse_model=CartesianDiscreteModel((0,1,0,1),(1,1)) |
| 140 | + dmodel=OctreeDistributedDiscreteModel(ranks,coarse_model,1) |
| 141 | + test_refine_and_coarsen_at_once(ranks,dmodel,order) |
| 142 | + rdmodel=dmodel |
| 143 | + for i=1:5 |
| 144 | + rdmodel=test_transfer_ops_and_redistribute(ranks,rdmodel,order) |
| 145 | + end |
| 146 | + end |
| 147 | + |
| 148 | + function test_3d(ranks,order) |
| 149 | + coarse_model=CartesianDiscreteModel((0,1,0,1,0,1),(2,2,2)) |
| 150 | + dmodel=OctreeDistributedDiscreteModel(ranks,coarse_model,0) |
| 151 | + test_refine_and_coarsen_at_once(ranks,dmodel,order) |
| 152 | + rdmodel=dmodel |
| 153 | + for i=1:5 |
| 154 | + rdmodel=test_transfer_ops_and_redistribute(ranks,rdmodel,order) |
| 155 | + end |
| 156 | + end |
| 157 | + |
| 158 | + u_ex_2D((x,y)) = 2*VectorValue(-y,x) |
| 159 | + f_ex_2D(x) = u_ex_2D(x) |
| 160 | + u_ex_3D((x,y,z)) = 2*VectorValue(-y,x,0.) - VectorValue(0.,-z,y) |
| 161 | + f_ex_3D(x) = u_ex_3D(x) |
| 162 | + |
| 163 | + function get_analytical_functions(Dc) |
| 164 | + if Dc==2 |
| 165 | + return u_ex_2D, f_ex_2D |
| 166 | + else |
| 167 | + @assert Dc==3 |
| 168 | + return u_ex_3D, f_ex_3D |
| 169 | + end |
| 170 | + end |
| 171 | + |
| 172 | + include("CoarseDiscreteModelsTools.jl") |
| 173 | + |
| 174 | + function solve_maxwell(model::GridapDistributed.DistributedDiscreteModel{Dc},order) where {Dc} |
| 175 | + u_ex, f_ex=get_analytical_functions(Dc) |
| 176 | + |
| 177 | + V = FESpace(model, |
| 178 | + ReferenceFE(nedelec,order), |
| 179 | + conformity=:Hcurl, |
| 180 | + dirichlet_tags="boundary") |
| 181 | + |
| 182 | + U = TrialFESpace(V,u_ex) |
| 183 | + |
| 184 | + trian = Triangulation(model) |
| 185 | + degree = 2*(order+1) |
| 186 | + dΩ = Measure(trian,degree) |
| 187 | + |
| 188 | + a(u,v) = ∫( (∇×u)⋅(∇×v) + u⋅v )dΩ |
| 189 | + l(v) = ∫(f_ex⋅v)dΩ |
| 190 | + |
| 191 | + op = AffineFEOperator(a,l,U,V) |
| 192 | + if (num_free_dofs(U)==0) |
| 193 | + # UMFPACK cannot handle empty linear systems |
| 194 | + uh = zero(U) |
| 195 | + else |
| 196 | + uh = solve(op) |
| 197 | + end |
| 198 | + |
| 199 | + # uh_ex=interpolate(u_ex_3D,U) |
| 200 | + # map(local_views(get_free_dof_values(uh_ex)), local_views(op.op.matrix), local_views(op.op.vector)) do U_ex, A, b |
| 201 | + # r_ex = A*U_ex - b |
| 202 | + # println(r_ex) |
| 203 | + # end |
| 204 | + uh,U |
| 205 | + end |
| 206 | + |
| 207 | + function check_error_maxwell(model::GridapDistributed.DistributedDiscreteModel{Dc},order,uh) where {Dc} |
| 208 | + trian = Triangulation(model) |
| 209 | + degree = 2*(order+1) |
| 210 | + dΩ = Measure(trian,degree) |
| 211 | + |
| 212 | + u_ex, f_ex = get_analytical_functions(Dc) |
| 213 | + |
| 214 | + eu = u_ex - uh |
| 215 | + |
| 216 | + l2(v) = sqrt(sum(∫(v⋅v)*dΩ)) |
| 217 | + hcurl(v) = sqrt(sum(∫(v⋅v + (∇×v)⋅(∇×v))*dΩ)) |
| 218 | + |
| 219 | + eu_l2 = l2(eu) |
| 220 | + eu_hcurl = hcurl(eu) |
| 221 | + |
| 222 | + tol = 1.0e-6 |
| 223 | + @test eu_l2 < tol |
| 224 | + @test eu_hcurl < tol |
| 225 | + end |
| 226 | + |
| 227 | + function run(distribute) |
| 228 | + # debug_logger = ConsoleLogger(stderr, Logging.Debug) |
| 229 | + # global_logger(debug_logger); # Enable the debug logger globally |
| 230 | + np = MPI.Comm_size(MPI.COMM_WORLD) |
| 231 | + ranks = distribute(LinearIndices((np,))) |
| 232 | + |
| 233 | + for order=0:2 |
| 234 | + test_2d(ranks,order) |
| 235 | + end |
| 236 | + |
| 237 | + for order=0:2 |
| 238 | + test_3d(ranks,order) |
| 239 | + end |
| 240 | + end |
| 241 | +end |
0 commit comments