Fixes

src/Solvers/gmg.jl

@@ -39,7 +39,7 @@ function gmg_solver(::Val{(:u,:j)},params)
         r = r + a_thin_wall_j_j(dj,v_j,p...)
       end
       if abs(ζ) > eps(typeof(ζ))
-        r = r + a_al_u_u(du,v_u,ζ,Πp,dΩf) + a_al_j_j(dj,v_j,ζ,dΩf)
+        r = r + a_al_sf(du,v_u,ζ,Πp,dΩf) + a_al_sf(dj,v_j,ζ,dΩf)
       end
       return r
     end
@@ -59,7 +59,7 @@ function gmg_solver(::Val{(:u,:j)},params)
         r = r + a_thin_wall_j_j(dj,v_j,p...)
       end
       if abs(ζ) > eps(typeof(ζ))
-        r = r + a_al_u_u(du,v_u,ζ,Πp,dΩf) + a_al_j_j(dj,v_j,ζ,dΩf)
+        r = r + a_al_sf(du,v_u,ζ,Πp,dΩf) + a_al_sf(dj,v_j,ζ,dΩf)
       end
       return r
     end
@@ -76,7 +76,7 @@ function gmg_solver(::Val{(:u,:j)},params)
         dc_mhd_u_u(u,du,v_u,α,dΩf)
       end
       if abs(ζ) > eps(typeof(ζ))
-        r = r + a_al_u_u(du,v_u,ζ,Πp,dΩf)
+        r = r + a_al_sf(du,v_u,ζ,Πp,dΩf)
       end
       return r
     end
@@ -87,14 +87,14 @@ function gmg_solver(::Val{(:u,:j)},params)
     _, _, α, β, γ, σf, f, B, ζ = retrieve_fluid_params(model,params,k)
     Ωf = Triangulation(model)
     dΩf = Measure(Ωf,2*k)
-    l((u,j),(v_u,v_j)) = a_al_u_u(u,v_u,ζ,Πp,dΩf) + a_al_j_j(j,v_j,ζ,dΩf)
+    l((u,j),(v_u,v_j)) = a_al_sf(u,v_u,ζ,Πp,dΩf) + a_al_sf(j,v_j,ζ,dΩf)
     return l
   end
   function rhs_u(model)
     _, _, α, β, γ, σf, f, B, ζ = retrieve_fluid_params(model,params,k)
     Ωf = Triangulation(model)
     dΩf = Measure(Ωf,2*k)
-    l(u,v_u) = a_al_u_u(u,v_u,ζ,Πp,dΩf)
+    l(u,v_u) = a_al_sf(u,v_u,ζ,Πp,dΩf)
     return l
   end
 

src/weakforms.jl

@@ -25,8 +25,7 @@ function _weak_form(params,k)
   bcs_params = retrieve_bcs_params(params,k)
   params_φ, params_thin_wall, params_f, params_B, params_Λ = bcs_params
 
-  reffe_p = params[:fespaces][:reffe_p]
-  Πp = MultilevelTools.LocalProjectionMap(divergence,reffe_p,2*k)
+  Πp = local_projection_operator(params,k)
 
   function a(x,dy)
     r = a_mhd(x,dy,β,γ,B,σf,dΩf)
@@ -105,8 +104,7 @@ function _ode_weak_form(params,k)
   bcs_params = retrieve_bcs_params(params,k)
   params_φ, params_thin_wall, params_f, params_B, params_Λ = bcs_params
 
-  reffe_p = params[:fespaces][:reffe_p]
-  Πp = MultilevelTools.LocalProjectionMap(divergence,reffe_p,2*k)
+  Πp = local_projection_operator(params,k)
 
   m(t,x,dy) = m_u(x,dy,dΩf)
 
@@ -330,10 +328,25 @@ p_dc_mhd_u_u(u,du,v_u,α,dΩ) = ∫( α*v_u⋅( conv∘(u,∇(du)) ) ) * dΩ
 function a_al(x,dy,ζ,Πp,dΩf,dΩ)
   u, p, j, φ = x
   v_u, v_p, v_j, v_φ = dy
-  a_al_u_u(u,v_u,ζ,Πp,dΩf) + a_al_j_j(j,v_j,ζ,dΩ)
+  if isnothing(Πp)
+    a_al_sf(u,v_u,ζ,dΩf) + a_al_sf(j,v_j,ζ,dΩ)
+  else
+    a_al_sf(u,v_u,ζ,Πp,dΩf) + a_al_sf(j,v_j,ζ,dΩ)
+  end
+end
+a_al_sf(x,y,ζ,Π,dΩ) = ∫(ζ*Π(x)*(∇⋅y))*dΩ
+a_al_sf(x,y,ζ,dΩ) = ∫(ζ*(∇⋅x)*(∇⋅y))*dΩ
+
+function local_projection_operator(params,k)
+  fluid_disc = params[:fespaces][:fluid_disc]
+  if fluid_disc ∈ (:SV,)
+    return nothing
+  end
+
+  reffe_p = params[:fespaces][:reffe_p]
+  Πp = MultilevelTools.LocalProjectionMap(divergence,reffe_p,2*k)
+  return Πp
 end
-a_al_u_u(u,v_u,ζ,Πp,dΩ) = ∫( ζ*Πp(u)*(∇⋅v_u) )*dΩ
-a_al_j_j(j,v_j,ζ,dΩ) = ∫( ζ*(∇⋅j)*(∇⋅v_j) )*dΩ
 
 # Solid equations