!-------------------------------------------------------------------------- !---------------------------- caldyn_slow_NH ---------------------------------- !$OMP DO SCHEDULE(STATIC) DO ij = 1, primal_num kdown = 1 kup = 1 l=1 w_il(l,ij) = 2.*W(l,ij)/rhodz(kup,ij) !DIR$ SIMD DO l = 2, llm kdown = l-1 kup = l w_il(l,ij) = 2.*W(l,ij)/(rhodz(kdown,ij)+rhodz(kup,ij)) END DO kdown = llm kup = llm l=llm+1 w_il(l,ij) = 2.*W(l,ij)/rhodz(kdown,ij) END DO !$OMP END DO !$OMP DO SCHEDULE(STATIC) DO edge = 1, edge_num ij_left = left(edge) ij_right = right(edge) kdown = 1 kup = 1 l=1 ! compute DePhi, v_el, G_el, F_el ! v_el, W2_el and therefore G_el incorporate metric factor le_de ! while DePhil, W_el and F_el do not W_el = .5*( W(l,ij_right)+W(l,ij_left) ) DePhil(l,edge) = 1.*(Phi(l,ij_right)-Phi(l,ij_left)) F_el(l,edge) = DePhil(l,edge)*W_el W2_el = .5*le_de(edge) * ( W(l,ij_left)*w_il(l,ij_left) + W(l,ij_right)*w_il(l,ij_right) ) v_el(l,edge) = .5*le_de(edge)*(u(kup,edge)+u(kdown,edge)) ! checked G_el(l,edge) = v_el(l,edge)*W_el - DePhil(l,edge)*W2_el !DIR$ SIMD DO l = 2, llm kdown = l-1 kup = l ! compute DePhi, v_el, G_el, F_el ! v_el, W2_el and therefore G_el incorporate metric factor le_de ! while DePhil, W_el and F_el do not W_el = .5*( W(l,ij_right)+W(l,ij_left) ) DePhil(l,edge) = 1.*(Phi(l,ij_right)-Phi(l,ij_left)) F_el(l,edge) = DePhil(l,edge)*W_el W2_el = .5*le_de(edge) * ( W(l,ij_left)*w_il(l,ij_left) + W(l,ij_right)*w_il(l,ij_right) ) v_el(l,edge) = .5*le_de(edge)*(u(kup,edge)+u(kdown,edge)) ! checked G_el(l,edge) = v_el(l,edge)*W_el - DePhil(l,edge)*W2_el END DO kdown = llm kup = llm l=llm+1 ! compute DePhi, v_el, G_el, F_el ! v_el, W2_el and therefore G_el incorporate metric factor le_de ! while DePhil, W_el and F_el do not W_el = .5*( W(l,ij_right)+W(l,ij_left) ) DePhil(l,edge) = 1.*(Phi(l,ij_right)-Phi(l,ij_left)) F_el(l,edge) = DePhil(l,edge)*W_el W2_el = .5*le_de(edge) * ( W(l,ij_left)*w_il(l,ij_left) + W(l,ij_right)*w_il(l,ij_right) ) v_el(l,edge) = .5*le_de(edge)*(u(kup,edge)+u(kdown,edge)) ! checked G_el(l,edge) = v_el(l,edge)*W_el - DePhil(l,edge)*W2_el END DO !$OMP END DO ! compute GradPhi2, dPhi, dW !$OMP DO SCHEDULE(STATIC) DO ij = 1, primal_num ! this VLOOP iterates over primal cell edges SELECT CASE(primal_deg(ij)) CASE(4) edge1 = primal_edge(1,ij) edge2 = primal_edge(2,ij) edge3 = primal_edge(3,ij) edge4 = primal_edge(4,ij) le_de1 = le_de(edge1) le_de2 = le_de(edge2) le_de3 = le_de(edge3) le_de4 = le_de(edge4) sign1 = primal_ne(1,ij) sign2 = primal_ne(2,ij) sign3 = primal_ne(3,ij) sign4 = primal_ne(4,ij) !DIR$ SIMD DO l = 1, llm+1 gPhi2=0. dP=0. divG=0 gPhi2 = gPhi2 + le_de1*DePhil(l,edge1)**2 dP = dP + le_de1*DePhil(l,edge1)*v_el(l,edge1) divG = divG + sign1*G_el(l,edge1) ! -div(G_el), G_el already has le_de gPhi2 = gPhi2 + le_de2*DePhil(l,edge2)**2 dP = dP + le_de2*DePhil(l,edge2)*v_el(l,edge2) divG = divG + sign2*G_el(l,edge2) ! -div(G_el), G_el already has le_de gPhi2 = gPhi2 + le_de3*DePhil(l,edge3)**2 dP = dP + le_de3*DePhil(l,edge3)*v_el(l,edge3) divG = divG + sign3*G_el(l,edge3) ! -div(G_el), G_el already has le_de gPhi2 = gPhi2 + le_de4*DePhil(l,edge4)**2 dP = dP + le_de4*DePhil(l,edge4)*v_el(l,edge4) divG = divG + sign4*G_el(l,edge4) ! -div(G_el), G_el already has le_de gradPhi2(l,ij) = 1./(2.*Ai(ij)) * gPhi2 dPhi(l,ij) = gradPhi2(l,ij)*w_il(l,ij) - 1./(2.*Ai(ij))*dP dW(l,ij) = (-1./Ai(ij))*divG END DO CASE(5) edge1 = primal_edge(1,ij) edge2 = primal_edge(2,ij) edge3 = primal_edge(3,ij) edge4 = primal_edge(4,ij) edge5 = primal_edge(5,ij) le_de1 = le_de(edge1) le_de2 = le_de(edge2) le_de3 = le_de(edge3) le_de4 = le_de(edge4) le_de5 = le_de(edge5) sign1 = primal_ne(1,ij) sign2 = primal_ne(2,ij) sign3 = primal_ne(3,ij) sign4 = primal_ne(4,ij) sign5 = primal_ne(5,ij) !DIR$ SIMD DO l = 1, llm+1 gPhi2=0. dP=0. divG=0 gPhi2 = gPhi2 + le_de1*DePhil(l,edge1)**2 dP = dP + le_de1*DePhil(l,edge1)*v_el(l,edge1) divG = divG + sign1*G_el(l,edge1) ! -div(G_el), G_el already has le_de gPhi2 = gPhi2 + le_de2*DePhil(l,edge2)**2 dP = dP + le_de2*DePhil(l,edge2)*v_el(l,edge2) divG = divG + sign2*G_el(l,edge2) ! -div(G_el), G_el already has le_de gPhi2 = gPhi2 + le_de3*DePhil(l,edge3)**2 dP = dP + le_de3*DePhil(l,edge3)*v_el(l,edge3) divG = divG + sign3*G_el(l,edge3) ! -div(G_el), G_el already has le_de gPhi2 = gPhi2 + le_de4*DePhil(l,edge4)**2 dP = dP + le_de4*DePhil(l,edge4)*v_el(l,edge4) divG = divG + sign4*G_el(l,edge4) ! -div(G_el), G_el already has le_de gPhi2 = gPhi2 + le_de5*DePhil(l,edge5)**2 dP = dP + le_de5*DePhil(l,edge5)*v_el(l,edge5) divG = divG + sign5*G_el(l,edge5) ! -div(G_el), G_el already has le_de gradPhi2(l,ij) = 1./(2.*Ai(ij)) * gPhi2 dPhi(l,ij) = gradPhi2(l,ij)*w_il(l,ij) - 1./(2.*Ai(ij))*dP dW(l,ij) = (-1./Ai(ij))*divG END DO CASE(6) edge1 = primal_edge(1,ij) edge2 = primal_edge(2,ij) edge3 = primal_edge(3,ij) edge4 = primal_edge(4,ij) edge5 = primal_edge(5,ij) edge6 = primal_edge(6,ij) le_de1 = le_de(edge1) le_de2 = le_de(edge2) le_de3 = le_de(edge3) le_de4 = le_de(edge4) le_de5 = le_de(edge5) le_de6 = le_de(edge6) sign1 = primal_ne(1,ij) sign2 = primal_ne(2,ij) sign3 = primal_ne(3,ij) sign4 = primal_ne(4,ij) sign5 = primal_ne(5,ij) sign6 = primal_ne(6,ij) !DIR$ SIMD DO l = 1, llm+1 gPhi2=0. dP=0. divG=0 gPhi2 = gPhi2 + le_de1*DePhil(l,edge1)**2 dP = dP + le_de1*DePhil(l,edge1)*v_el(l,edge1) divG = divG + sign1*G_el(l,edge1) ! -div(G_el), G_el already has le_de gPhi2 = gPhi2 + le_de2*DePhil(l,edge2)**2 dP = dP + le_de2*DePhil(l,edge2)*v_el(l,edge2) divG = divG + sign2*G_el(l,edge2) ! -div(G_el), G_el already has le_de gPhi2 = gPhi2 + le_de3*DePhil(l,edge3)**2 dP = dP + le_de3*DePhil(l,edge3)*v_el(l,edge3) divG = divG + sign3*G_el(l,edge3) ! -div(G_el), G_el already has le_de gPhi2 = gPhi2 + le_de4*DePhil(l,edge4)**2 dP = dP + le_de4*DePhil(l,edge4)*v_el(l,edge4) divG = divG + sign4*G_el(l,edge4) ! -div(G_el), G_el already has le_de gPhi2 = gPhi2 + le_de5*DePhil(l,edge5)**2 dP = dP + le_de5*DePhil(l,edge5)*v_el(l,edge5) divG = divG + sign5*G_el(l,edge5) ! -div(G_el), G_el already has le_de gPhi2 = gPhi2 + le_de6*DePhil(l,edge6)**2 dP = dP + le_de6*DePhil(l,edge6)*v_el(l,edge6) divG = divG + sign6*G_el(l,edge6) ! -div(G_el), G_el already has le_de gradPhi2(l,ij) = 1./(2.*Ai(ij)) * gPhi2 dPhi(l,ij) = gradPhi2(l,ij)*w_il(l,ij) - 1./(2.*Ai(ij))*dP dW(l,ij) = (-1./Ai(ij))*divG END DO CASE DEFAULT !DIR$ SIMD DO l = 1, llm+1 gPhi2=0. dP=0. divG=0 DO iedge = 1, primal_deg(ij) edge = primal_edge(iedge,ij) gPhi2 = gPhi2 + le_de(edge)*DePhil(l,edge)**2 dP = dP + le_de(edge)*DePhil(l,edge)*v_el(l,edge) divG = divG + primal_ne(iedge,ij)*G_el(l,edge) ! -div(G_el), G_el already has le_de END DO gradPhi2(l,ij) = 1./(2.*Ai(ij)) * gPhi2 dPhi(l,ij) = gradPhi2(l,ij)*w_il(l,ij) - 1./(2.*Ai(ij))*dP dW(l,ij) = (-1./Ai(ij))*divG END DO END SELECT END DO !$OMP END DO ! We need a barrier here because we compute gradPhi2, F_el and w_il above and do a vertical average below !$OMP BARRIER ! Compute berni at scalar points !$OMP DO SCHEDULE(STATIC) DO ij = 1, primal_num ! this VLOOP iterates over primal cell edges SELECT CASE(primal_deg(ij)) CASE(4) edge1 = primal_edge(1,ij) edge2 = primal_edge(2,ij) edge3 = primal_edge(3,ij) edge4 = primal_edge(4,ij) le_de1 = le_de(edge1) le_de2 = le_de(edge2) le_de3 = le_de(edge3) le_de4 = le_de(edge4) !DIR$ SIMD DO l = 1, llm u2=0. u2 = u2 + le_de1*u(l,edge1)**2 u2 = u2 + le_de2*u(l,edge2)**2 u2 = u2 + le_de3*u(l,edge3)**2 u2 = u2 + le_de4*u(l,edge4)**2 berni(l,ij) = 1./(4.*Ai(ij)) * u2 - .25*( gradPhi2(l,ij)*w_il(l,ij)**2 + gradPhi2(l+1,ij)*w_il(l+1,ij)**2 ) END DO CASE(5) edge1 = primal_edge(1,ij) edge2 = primal_edge(2,ij) edge3 = primal_edge(3,ij) edge4 = primal_edge(4,ij) edge5 = primal_edge(5,ij) le_de1 = le_de(edge1) le_de2 = le_de(edge2) le_de3 = le_de(edge3) le_de4 = le_de(edge4) le_de5 = le_de(edge5) !DIR$ SIMD DO l = 1, llm u2=0. u2 = u2 + le_de1*u(l,edge1)**2 u2 = u2 + le_de2*u(l,edge2)**2 u2 = u2 + le_de3*u(l,edge3)**2 u2 = u2 + le_de4*u(l,edge4)**2 u2 = u2 + le_de5*u(l,edge5)**2 berni(l,ij) = 1./(4.*Ai(ij)) * u2 - .25*( gradPhi2(l,ij)*w_il(l,ij)**2 + gradPhi2(l+1,ij)*w_il(l+1,ij)**2 ) END DO CASE(6) edge1 = primal_edge(1,ij) edge2 = primal_edge(2,ij) edge3 = primal_edge(3,ij) edge4 = primal_edge(4,ij) edge5 = primal_edge(5,ij) edge6 = primal_edge(6,ij) le_de1 = le_de(edge1) le_de2 = le_de(edge2) le_de3 = le_de(edge3) le_de4 = le_de(edge4) le_de5 = le_de(edge5) le_de6 = le_de(edge6) !DIR$ SIMD DO l = 1, llm u2=0. u2 = u2 + le_de1*u(l,edge1)**2 u2 = u2 + le_de2*u(l,edge2)**2 u2 = u2 + le_de3*u(l,edge3)**2 u2 = u2 + le_de4*u(l,edge4)**2 u2 = u2 + le_de5*u(l,edge5)**2 u2 = u2 + le_de6*u(l,edge6)**2 berni(l,ij) = 1./(4.*Ai(ij)) * u2 - .25*( gradPhi2(l,ij)*w_il(l,ij)**2 + gradPhi2(l+1,ij)*w_il(l+1,ij)**2 ) END DO CASE DEFAULT !DIR$ SIMD DO l = 1, llm u2=0. DO iedge = 1, primal_deg(ij) edge = primal_edge(iedge,ij) u2 = u2 + le_de(edge)*u(l,edge)**2 END DO berni(l,ij) = 1./(4.*Ai(ij)) * u2 - .25*( gradPhi2(l,ij)*w_il(l,ij)**2 + gradPhi2(l+1,ij)*w_il(l+1,ij)**2 ) END DO END SELECT END DO !$OMP END DO !$OMP DO SCHEDULE(STATIC) DO edge = 1, edge_num ij_left = left(edge) ij_right = right(edge) !DIR$ SIMD DO l = 1, llm ! Compute mass flux and grad(berni) uu = .5*(rhodz(l,ij_left)+rhodz(l,ij_right))*u(l,edge) - .5*( F_el(l,edge)+F_el(l+1,edge) ) hflux(l,edge) = le_de(edge)*uu du(l,edge) = 1.*(berni(l,ij_left)-berni(l,ij_right)) END DO END DO !$OMP END DO !---------------------------- caldyn_slow_NH ---------------------------------- !--------------------------------------------------------------------------