Merge lp:~s-parkinson11/fluidity/SU_SUPG_stress_form_fix into lp:fluidity

[Samuel Parkinson]

[Stephan Kramer (s-kramer)]

Looks good to me. The only issue that one could raise is that this only deals with isotropic viscosity in full/partial stress form - but in my entirely uninformed opinion we don't deal with anisotropic full stress form very well in general. So afaic this is good to go when tests have gone green. Oh, one nitpicky thing: I much prefer ==, >=, <, etc. over .eq., .geq. and .lt.

review: Approve

=== modified file 'assemble/Momentum_CG.F90'
--- assemble/Momentum_CG.F90	2012-05-23 14:01:16 +0000
+++ assemble/Momentum_CG.F90	2012-08-17 13:45:11 +0000
@@ -1184,7 +1184,7 @@
       real, dimension(u%dim, ele_loc(u, ele)) :: oldu_val
       type(element_type), pointer :: u_shape, p_shape
       real, dimension(ele_ngi(u, ele)) :: detwei, detwei_old, detwei_new
-      real, dimension(u%dim, u%dim, ele_ngi(u,ele)) :: J_mat
+      real, dimension(u%dim, u%dim, ele_ngi(u,ele)) :: J_mat, diff_q
       real, dimension(ele_loc(u, ele), ele_ngi(u, ele), u%dim) :: du_t
       real, dimension(ele_loc(u, ele), ele_ngi(u, ele), u%dim) :: dug_t
       real, dimension(ele_loc(p, ele), ele_ngi(p, ele), u%dim) :: dp_t
@@ -1197,7 +1197,7 @@
       real, dimension(u%dim, ele_loc(u, ele)) :: big_m_diag_addto, rhs_addto
       real, dimension(u%dim, u%dim, ele_loc(u, ele), ele_loc(u, ele)) :: big_m_tensor_addto
       logical, dimension(u%dim, u%dim) :: block_mask ! control whether the off diagonal entries are used
-      integer :: dim
+      integer :: dim, i, j
       type(element_type) :: test_function
 
       if(move_mesh) then
@@ -1272,11 +1272,22 @@
             relu_gi = relu_gi - ele_val_at_quad(ug, ele)
           end if
           if(have_viscosity) then
-          call supg_test_function(supg_shape, u_shape, du_t, relu_gi, j_mat, diff_q = ele_val_at_quad(viscosity, ele), &
-            & nu_bar_scheme = nu_bar_scheme, nu_bar_scale = nu_bar_scale)
+             diff_q = ele_val_at_quad(viscosity, ele)
+
+             ! for full and partial stress form we need to set the off diagonal terms of the viscosity tensor to zero
+             ! to be able to invert it when calculating nu_bar
+             do i=1,size(diff_q,1)
+                do j=1,size(diff_q,2)
+                   if(i.eq.j) cycle
+                   diff_q(i,j,:) = 0.0
+                end do
+             end do
+
+             call supg_test_function(supg_shape, u_shape, du_t, relu_gi, j_mat, diff_q = ele_val_at_quad(viscosity, ele), &
+                  & nu_bar_scheme = nu_bar_scheme, nu_bar_scale = nu_bar_scale)
           else
-          call supg_test_function(supg_shape, u_shape, du_t, relu_gi, j_mat, &
-            & nu_bar_scheme = nu_bar_scheme, nu_bar_scale = nu_bar_scale)
+             call supg_test_function(supg_shape, u_shape, du_t, relu_gi, j_mat, &
+                  & nu_bar_scheme = nu_bar_scheme, nu_bar_scale = nu_bar_scale)
           end if
           test_function = supg_shape
         case default
@@ -1531,11 +1542,11 @@
       real, dimension(ele_loc(u, ele), ele_ngi(u, ele), u%dim), intent(in) :: dug_t
       real, dimension(:, :, :), intent(in) :: dnvfrac_t
       real, dimension(ele_ngi(u, ele)), intent(in) :: detwei
-      real, dimension(u%dim, u%dim, ele_ngi(u,ele)) :: J_mat
+      real, dimension(u%dim, u%dim, ele_ngi(u,ele)) :: J_mat, diff_q
       real, dimension(u%dim, u%dim, ele_loc(u, ele), ele_loc(u, ele)), intent(inout) :: big_m_tensor_addto
       real, dimension(u%dim, ele_loc(u, ele)), intent(inout) :: rhs_addto
     
-      integer :: dim, i
+      integer :: dim, i, j
       real, dimension(ele_ngi(u, ele)) :: density_gi, div_relu_gi
       real, dimension(ele_ngi(u, ele)) :: nvfrac_gi, relu_dot_grad_nvfrac_gi
       real, dimension(u%dim, ele_ngi(u, ele)) :: grad_nvfrac_gi
@@ -1603,15 +1614,26 @@
       
       select case(stabilisation_scheme)
       case(STABILISATION_STREAMLINE_UPWIND)
-        if(have_viscosity) then
-          advection_mat = advection_mat + &
-            & element_upwind_stabilisation(u_shape, du_t, relu_gi, J_mat, detwei, &
-            & diff_q = ele_val_at_quad(viscosity, ele), nu_bar_scheme = nu_bar_scheme, nu_bar_scale = nu_bar_scale)
-        else
-           advection_mat = advection_mat + &
-            & element_upwind_stabilisation(u_shape, du_t, relu_gi, J_mat, detwei, &
-            & nu_bar_scheme = nu_bar_scheme, nu_bar_scale = nu_bar_scale)
-        end if
+         if(have_viscosity) then
+            diff_q = ele_val_at_quad(viscosity, ele)
+
+            ! for full and partial stress form we need to set the off diagonal terms of the viscosity tensor to zero
+            ! to be able to invert it when calculating nu_bar
+            do i=1,size(diff_q,1)
+               do j=1,size(diff_q,2)
+                  if(i.eq.j) cycle
+                  diff_q(i,j,:) = 0.0
+               end do
+            end do
+
+            advection_mat = advection_mat + &
+                 & element_upwind_stabilisation(u_shape, du_t, relu_gi, J_mat, detwei, &
+                 & diff_q, nu_bar_scheme = nu_bar_scheme, nu_bar_scale = nu_bar_scale)
+         else
+            advection_mat = advection_mat + &
+                 & element_upwind_stabilisation(u_shape, du_t, relu_gi, J_mat, detwei, &
+                 & nu_bar_scheme = nu_bar_scheme, nu_bar_scale = nu_bar_scale)
+         end if
       end select
       
       do dim = 1, u%dim

=== modified file 'femtools/Diagnostic_Fields.F90'
--- femtools/Diagnostic_Fields.F90	2012-02-18 15:45:48 +0000
+++ femtools/Diagnostic_Fields.F90	2012-08-17 13:45:11 +0000
@@ -477,7 +477,7 @@
     type(scalar_field), intent(inout) :: grn
 
     type(vector_field), pointer :: U, X
-    integer :: ele, gi, stat
+    integer :: ele, gi, stat, a, b
     ! Transformed quadrature weights.
     real, dimension(ele_ngi(GRN, 1)) :: detwei
     ! Inverse of the local coordinate change matrix.
@@ -525,6 +525,24 @@
        ! The matmul is as given by dham
        GRN_q=ele_val_at_quad(U, ele)
        vis_q=ele_val_at_quad(viscosity, ele)
+
+       ! for full and partial stress form we need to set the off diagonal terms of the viscosity tensor to zero
+       ! to be able to invert it 
+       if (have_option(trim(U%option_path)//&
+            &"/prognostic/spatial_discretisation/continuous_galerkin"//&
+            &"/stress_terms/stress_form") .or. &
+            have_option(trim(U%option_path)//&
+            &"/prognostic/spatial_discretisation/continuous_galerkin"//&
+            &"/stress_terms/partial_stress_form")) then
+
+          do a=1,size(vis_q,1)
+             do b=1,size(vis_q,2)
+                if(a.eq.b) cycle
+                vis_q(a,b,:) = 0.0
+             end do
+          end do
+       end if
+
        do gi=1, size(detwei)
           GRN_q(:,gi)=matmul(GRN_q(:,gi), J(:,:,gi))
           GRN_q(:,gi)=matmul(inverse(vis_q(:,:,gi)), GRN_q(:,gi))