gvec_post.F90

!!matvec with matmul !#define __MATVEC_N(y,Mat,Vec) y=MATMUL(Mat,Vec) !#define __MATVEC_T(y,Mat,Vec) y=MATMUL(Vec,Mat) !#define __PMATVEC_N(fy,y,Mat,Vec) y=fyy+MATMUL(Mat,Vec) !#define __PMATVEC_T(fy,y,Mat,Vec) y=fyy+MATMUL(Vec,Mat) !#define __AMATVEC_N(y,fMat,Mat,Vec) y=fMatMATMUL(Mat,Vec) !#define __AMATVEC_T(y,fMat,Mat,Vec) y=fMatMATMUL(Vec,Mat) !#define __PAMATVEC_N(fy,y,fMat,Mat,Vec) y=fyy+fMatMATMUL(Mat,Vec) !#define __PAMATVEC_T(fy,y,fMat,Mat,Vec) y=fyy+fMatMATMUL(Vec,Mat)

!!#define __GENERICMATVEC(NT,fy,y,fMat,Mat,Vec) CALL DGEMV(NT,SIZE(Mat,1),SIZE(Mat,2),fMat,Mat,SIZE(Mat,1),Vec,1,fy,y,1)

!!matmat with matmul !#define __MATMAT_NN(Y,A,B) Y=MATMUL(A,B) !#define __MATMAT_TN(Y,A,B) Y=MATMUL(TRANSPOSE(A),B) !#define __MATMAT_NT(Y,A,B) Y=MATMUL(A,TRANSPOSE(B)) !#define __MATMAT_TT(Y,A,B) Y=TRANSPOSE(MATMUL(B,A))

!#define __PMATMAT_NN(fy,Y,A,B) Y=fyY+MATMUL(A,B) !#define __PMATMAT_TN(fy,Y,A,B) Y=fyY+MATMUL(TRANSPOSE(A),B) !#define __PMATMAT_NT(fy,Y,A,B) Y=fyY+MATMUL(A,TRANSPOSE(B)) !#define __PMATMAT_TT(fy,Y,A,B) Y=fyY+TRANSPOSE(MATMUL(B,A))

!#define __AMATMAT_NN(Y,fa,A,B) Y=faMATMUL(A,B) !#define __AMATMAT_TN(Y,fa,A,B) Y=faMATMUL(TRANSPOSE(A),B) !#define __AMATMAT_NT(Y,fa,A,B) Y=faMATMUL(A,TRANSPOSE(B)) !#define __AMATMAT_TT(Y,fa,A,B) Y=faTRANSPOSE(MATMUL(B,A))

!#define __PAMATMAT_NN(fy,Y,fa,A,B) Y=fyY+faMATMUL(A,B) !#define __PAMATMAT_TN(fy,Y,fa,A,B) Y=fyY+faMATMUL(TRANSPOSE(A),B) !#define __PAMATMAT_NT(fy,Y,fa,A,B) Y=fyY+faMATMUL(A,TRANSPOSE(B)) !#define __PAMATMAT_TT(fy,Y,fa,A,B) Y=fyY+faTRANSPOSE(MATMUL(B,A))

! GEMM does in general Y = fa A^?B^? + fy Y ! with structure: (m x n) = (m x k) (k x n) ! Y=A B : DGEMM('N','N',m,n,k,fa,Amat ,m, Bmat,k, fy,Y,m) ! Y=A^TB : DGEMM('T','N',m,n,k,fa,Amat ,k, Bmat,k, fy,Y,m) ! Y=A B^T : DGEMM('N','T',m,n,k,fa,Amat ,m, Bmat,n, fy,Y,m) ! Y=A^T*B^T : DGEMM('T','T',m,n,k,fa,Amat ,k, Bmat,n, fy,Y,m)

!#define __GENERICMATMAT_NN(fy,Y,fa,A,B) CALL DGEMM('N','N',SIZE(A,1),SIZE(B,2),SIZE(B,1),fa,A,SIZE(A,1),B,SIZE(B,1),fy,Y,SIZE(A,1)) !#define __GENERICMATMAT_TN(fy,Y,fa,A,B) CALL DGEMM('T','N',SIZE(A,2),SIZE(B,2),SIZE(B,1),fa,A,SIZE(B,1),B,SIZE(B,1),fy,Y,SIZE(A,2)) !#define __GENERICMATMAT_NT(fy,Y,fa,A,B) CALL DGEMM('N','T',SIZE(A,1),SIZE(B,1),SIZE(B,2),fa,A,SIZE(A,1),B,SIZE(B,1),fy,Y,SIZE(A,1)) !#define __GENERICMATMAT_TT(fy,Y,fa,A,B) CALL DGEMM('T','T',SIZE(A,2),SIZE(B,1),SIZE(B,2),fa,A,SIZE(B,2),B,SIZE(B,1),fy,Y,SIZE(A,2))

! SIMPLE INTERFACE FOR DGEMM, specifying nrows/ncols of mat A and nrows/ncols of mat B (for any transpose!) ! GEMM does in general Y = fa A^?B^? + fy Y ! with structure: (m x n) = (m x k) (k x n) ! Y=A B : DGEMM('N','N',m,n,k,fa,Amat ,m, Bmat,k, fy,Y,m) ! Y=A^TB : DGEMM('T','N',m,n,k,fa,Amat ,k, Bmat,k, fy,Y,m) ! Y=A B^T : DGEMM('N','T',m,n,k,fa,Amat ,m, Bmat,n, fy,Y,m) ! Y=A^T*B^T : DGEMM('T','T',m,n,k,fa,Amat ,k, Bmat,n, fy,Y,m)

Source Code

!===================================================================================================================================
! Copyright (c) 2025 GVEC Contributors, Max Planck Institute for Plasma Physics
! License: MIT
!===================================================================================================================================
#include "defines.FPP"


!===================================================================================================================================
!>
!!# **GVEC** Driver program
!!
!===================================================================================================================================
PROGRAM GVEC_POST
  USE MODgvec_MPI          ,ONLY : par_Init, par_Finalize,par_bcast
  USE MODgvec_Globals
  USE MODgvec_Analyze      ,ONLY: InitAnalyze,Analyze,FinalizeAnalyze
  USE MODgvec_Output       ,ONLY: InitOutput,FinalizeOutput
  USE MODgvec_Output_vars  ,ONLY: OutputLevel
  USE MODgvec_Restart      ,ONLY: InitRestart,FinalizeRestart
  USE MODgvec_Restart      ,ONLY: RestartFromState
  USE MODgvec_Output_Vars  ,ONLY: OutputLevel,ProjectName
  USE MODgvec_ReadState_Vars,ONLY: fileID_r,outputLevel_r
  USE MODgvec_MHD3D_visu   ,ONLY:WriteSFLoutfile
  USE MODgvec_MHD3D_EvalFunc , ONLY: InitProfilesGP,EvalEnergy,EvalForce
  USE MODgvec_ReadInTools  ,ONLY: FillStrings,GETLOGICAL,GETINT,IgnoredStrings
  USE MODgvec_MHD3D, ONLY:t_functional_mhd3d
!$ USE omp_lib
  IMPLICIT NONE
!-----------------------------------------------------------------------------------------------------------------------------------
!local variables
  INTEGER                 :: iArg,nArgs
  CHARACTER(LEN=255)      :: Parameterfile
  CHARACTER(LEN=255)      :: Statefile
  INTEGER                 :: JacCheck
  CLASS(t_functional_mhd3d),ALLOCATABLE :: functional
  REAL(wp)                :: StartTime,EndTime
!===================================================================================================================================
  CALL par_Init()
  __PERFINIT
  __PERFON('main')
  nArgs=COMMAND_ARGUMENT_COUNT()
  IF ((nArgs.LT.2))THEN
    ! Print out error message containing valid syntax
    STOP 'ERROR - Invalid syntax. Please use: gvec_post parameter.ini [Statefiles*] '
  END IF
  CALL GET_COMMAND_ARGUMENT(1,Parameterfile)

  CALL CPU_TIME(StartTime)
!$ StartTime=OMP_GET_WTIME()

  !header
  SWRITE(Unit_stdOut,'(132("="))')
  SWRITE(UNIT_stdOut,'(A)') "GVEC POST ! GVEC POST ! GVEC POST ! GVEC POST"
  SWRITE(Unit_stdOut,'(132("="))')
  !.only executes if compiled with OpenMP
!$ SWRITE(UNIT_stdOut,'(A,I6)')'   Number of OpenMP threads : ',OMP_GET_MAX_THREADS()
!$ SWRITE(Unit_stdOut,'(132("="))')
  !.only executes if compiled with MPI
# if MPI
  SWRITE(UNIT_stdOut,'(A,I6)')'   Number of MPI tasks : ',nRanks
  SWRITE(Unit_stdOut,fmt_sep)
  IF(nRanks.GT.1) CALL abort(__STAMP__,&
                   "GVEC post is compiled with MPI, but can only be called with 1 MPI rank." )
# endif
#include  "configuration-cmake.F90"
  SWRITE(Unit_stdOut,fmt_sep)

  CALL FillStrings(ParameterFile) !< readin parameterfile, done on MPI root + Bcast
  testdbg =.FALSE.
  testlevel=-1

  !initialization phase
  CALL InitRestart()
  CALL InitOutput()
  CALL InitAnalyze()

  ALLOCATE(t_functional_mhd3d :: functional)
  CALL functional%init()

  CALL IgnoredStrings()
  DO iArg=2,nArgs
    CALL GET_COMMAND_ARGUMENT(iArg,StateFile)
    SWRITE(Unit_stdOut,'(132("-"))')
    SWRITE(UNIT_stdOut,'(A,I4,A4,I4,A3,A)') 'Post-Analyze StateFile ',iArg-1,' of ',nArgs-1,' : ',TRIM(StateFile)
    SWRITE(Unit_stdOut,'(132("-"))')
    ProjectName='POST_'//TRIM(StateFile(1:INDEX(StateFile,'_State_')-1))
    ASSOCIATE(vars=>functional%minimizer%vars)
      CALL RestartFromState(StateFile,vars%dofs(0))
      outputLevel=outputLevel_r
      JacCheck=2
      !...check this: temporarily commented for gvec_post to run with MPI version...
      CALL InitProfilesGP() !evaluate profiles once at Gauss Points (on MPIroot + BCast)
      vars%dofs(0)%W_MHD3D=EvalEnergy(vars%dofs(0),.TRUE.,JacCheck)
      CALL EvalForce(vars%dofs(0),.FALSE.,JacCheck, vars%force(0))
      CALL Analyze(FileID_r, vars%dofs(0), vars%force(0))
      CALL writeSFLoutfile(vars%dofs(0),FileID_r)
    END ASSOCIATE
  END DO !iArg

  IF(ALLOCATED(functional))THEN
    CALL functional%free()
    DEALLOCATE(functional)
  END IF

  CALL FinalizeAnalyze()
  CALL FinalizeOutput()
  CALL FinalizeRestart()

  CALL CPU_TIME(EndTime)
!$ EndTime=OMP_GET_WTIME()
  WRITE(Unit_stdOut,'(132("="))')
  WRITE(Unit_stdOut,'(A,F8.2,A)') ' GVEC POST FINISHED ! [',EndTime-StartTime,' sec ]'
  WRITE(Unit_stdOut,'(132("="))')

  __PERFOFF('main')
  __PERFOUT('main')
  CALL par_Finalize()

END PROGRAM GVEC_POST

gvec_post.F90 Source File

Contents

Programs

Source Code

This file depends on

Source Code