out stability metrics to h5 file

PrincetonUniversity · Mar 8, 2024 · 3c08d8f · 3c08d8f
1 parent 91dc66a
commit 3c08d8f
Show file tree

Hide file tree

Showing 2 changed files with 42 additions and 7 deletions.
diff --git a/src/hesian.f90 b/src/hesian.f90
@@ -37,7 +37,9 @@ subroutine hesian( NGdof, position, Mvol, mn, LGdof )
                         Lhessian3Dallocated,denergydrr, denergydrz,denergydzr,denergydzz, &
 					            	LocalConstraint
 
-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
+  use sphdf5, only : write_stability
+
+  !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
 
   LOCALS
 
@@ -418,12 +420,12 @@ subroutine hesian( NGdof, position, Mvol, mn, LGdof )
 
   if( LHmatrix ) then
 
-   if( myid.eq.0 ) then ; cput = GETTIME ; write(ounit,'("hesian : ",f10.2," : LHmatrix="L2" ;")')cput-cpus, LHmatrix ;
-    open(munit, file="."//trim(ext)//".GF.ma", status="unknown", form="unformatted")
-    write(munit) NGdof
-    write(munit) ohessian(1:NGdof,1:NGdof)
-    close(munit)
-   endif
+    if( myid.eq.0 ) then ; cput = GETTIME ; write(ounit,'("hesian : ",f10.2," : LHmatrix="L2" ;")')cput-cpus, LHmatrix ;
+      open(munit, file="."//trim(ext)//".GF.ma", status="unknown", form="unformatted")
+      write(munit) NGdof
+      write(munit) ohessian(1:NGdof,1:NGdof)
+      close(munit)
+    endif
 
   endif
 
@@ -517,6 +519,10 @@ subroutine hesian( NGdof, position, Mvol, mn, LGdof )
     enddo
    enddo
 
+  !  do ii = 1, mn ; write(ounit,*) 'psifact', ii, psifactor(ii,Mvol-1)
+  ! enddo
+
+
 !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
 
    if( myid.eq.0 ) then ! screen output; 04 Dec 14;
@@ -585,6 +591,9 @@ subroutine hesian( NGdof, position, Mvol, mn, LGdof )
 
   endif ! end of if(                 ( LHevalues .or. LHevectors ) )
 
+  ! output hessian, eigenvalues, and eigenvectors to the .h5 file
+  WCALL( hesian, write_stability, (ohessian, evalr, evali, evecr, NGdof) )
+
 !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
 
   if( myid.eq.0 .and. Lperturbed.eq.1 ) then

diff --git a/src/sphdf5.f90 b/src/sphdf5.f90
@@ -986,6 +986,32 @@ end subroutine write_vector_potential
 !> \brief Write the final state of the equilibrium to the output file.
 !> \ingroup grp_output
 !>
+
+subroutine write_stability(ohessian, evalr, evali, evecr, NGdof)
+
+  LOCALS
+  integer, intent(in) :: NGdof
+  REAL, intent(in)    :: ohessian(:,:), evalr(:), evali(:), evecr(:,:)
+  integer(hid_t)      :: grpStability
+
+  BEGIN( sphdf5 )
+
+ if (myid.eq.0 .and. .not.skip_write) then
+
+  HDEFGRP( file_id, stability, grpStability)
+
+  HWRITERA( grpStability, NGdof, NGdof, ohessian, ohessian(1:NGdof,1:NGdof) )
+  HWRITERA( grpStability, NGdof, NGdof, evecr, evecr(1:NGdof,1:NGdof) )
+  HWRITERV( grpStability, NGdof, evalr, evalr(1:NGdof) )
+  HWRITERV( grpStability, NGdof, evali, evali(1:NGdof) )
+
+  HCLOSEGRP( grpStability )
+
+ endif ! myid.eq.0
+
+end subroutine write_stability
+
+
 subroutine hdfint
 
   use fileunits, only : ounit