Write variances instead of RMSs, to avoid floating point exceptions.

CaNS-World · Dec 19, 2023 · 99778bb · 99778bb
1 parent b78414b
commit 99778bb
Showing 1 changed file with 16 additions and 16 deletions.
diff --git a/src/output.f90 b/src/output.f90
@@ -430,14 +430,14 @@ subroutine out2d_duct(fname,ng,lo,hi,idir,l,dl,z_g,u,v,w) ! e.g. for a duct
       call MPI_ALLREDUCE(MPI_IN_PLACE,w2(1,1),ng(1)*ng(3),MPI_REAL_RP,MPI_SUM,MPI_COMM_WORLD,ierr)
       call MPI_ALLREDUCE(MPI_IN_PLACE,vw(1,1),ng(1)*ng(3),MPI_REAL_RP,MPI_SUM,MPI_COMM_WORLD,ierr)
       call MPI_ALLREDUCE(MPI_IN_PLACE,uv(1,1),ng(1)*ng(3),MPI_REAL_RP,MPI_SUM,MPI_COMM_WORLD,ierr)
-      um(:,:) =      um(:,:)*grid_area_ratio
-      vm(:,:) =      vm(:,:)*grid_area_ratio
-      wm(:,:) =      wm(:,:)*grid_area_ratio
-      u2(:,:) = sqrt(u2(:,:)*grid_area_ratio - um(:,:)**2)
-      v2(:,:) = sqrt(v2(:,:)*grid_area_ratio - vm(:,:)**2)
-      w2(:,:) = sqrt(w2(:,:)*grid_area_ratio - wm(:,:)**2)
-      vw(:,:) =      vw(:,:)*grid_area_ratio - vm(:,:)*wm(:,:)
-      uv(:,:) =      uv(:,:)*grid_area_ratio - um(:,:)*vm(:,:)
+      um(:,:) = um(:,:)*grid_area_ratio
+      vm(:,:) = vm(:,:)*grid_area_ratio
+      wm(:,:) = wm(:,:)*grid_area_ratio
+      u2(:,:) = u2(:,:)*grid_area_ratio - um(:,:)**2
+      v2(:,:) = v2(:,:)*grid_area_ratio - vm(:,:)**2
+      w2(:,:) = w2(:,:)*grid_area_ratio - wm(:,:)**2
+      vw(:,:) = vw(:,:)*grid_area_ratio - vm(:,:)*wm(:,:)
+      uv(:,:) = uv(:,:)*grid_area_ratio - um(:,:)*vm(:,:)
       if(myid == 0) then
         open(newunit=iunit,file=fname)
         do k=1,ng(3)
@@ -488,14 +488,14 @@ subroutine out2d_duct(fname,ng,lo,hi,idir,l,dl,z_g,u,v,w) ! e.g. for a duct
       call MPI_ALLREDUCE(MPI_IN_PLACE,w2(1,1),ng(1)*ng(3),MPI_REAL_RP,MPI_SUM,MPI_COMM_WORLD,ierr)
       call MPI_ALLREDUCE(MPI_IN_PLACE,uv(1,1),ng(1)*ng(3),MPI_REAL_RP,MPI_SUM,MPI_COMM_WORLD,ierr)
       call MPI_ALLREDUCE(MPI_IN_PLACE,uw(1,1),ng(1)*ng(3),MPI_REAL_RP,MPI_SUM,MPI_COMM_WORLD,ierr)
-      um(:,:) =      um(:,:)*grid_area_ratio
-      vm(:,:) =      vm(:,:)*grid_area_ratio
-      wm(:,:) =      wm(:,:)*grid_area_ratio
-      u2(:,:) = sqrt(u2(:,:)*grid_area_ratio - um(:,:)**2)
-      v2(:,:) = sqrt(v2(:,:)*grid_area_ratio - vm(:,:)**2)
-      w2(:,:) = sqrt(w2(:,:)*grid_area_ratio - wm(:,:)**2)
-      uv(:,:) =      uv(:,:)*grid_area_ratio - um(:,:)*vm(:,:)
-      uw(:,:) =      uw(:,:)*grid_area_ratio - um(:,:)*wm(:,:)
+      um(:,:) = um(:,:)*grid_area_ratio
+      vm(:,:) = vm(:,:)*grid_area_ratio
+      wm(:,:) = wm(:,:)*grid_area_ratio
+      u2(:,:) = u2(:,:)*grid_area_ratio - um(:,:)**2
+      v2(:,:) = v2(:,:)*grid_area_ratio - vm(:,:)**2
+      w2(:,:) = w2(:,:)*grid_area_ratio - wm(:,:)**2
+      uv(:,:) = uv(:,:)*grid_area_ratio - um(:,:)*vm(:,:)
+      uw(:,:) = uw(:,:)*grid_area_ratio - um(:,:)*wm(:,:)
       if(myid == 0) then
         open(newunit=iunit,file=fname)
         do k=1,ng(3)