! $Id$ ! ! This module provide a way for users to specify custom ! (i.e. not in the standard Pencil Code) physics, diagnostics etc. ! ! The module provides a set of standard hooks into the Pencil-Code and ! currently allows the following customizations: ! ! Description | Relevant function call ! --------------------------------------------------------------------------- ! Special variable registration | register_special ! (pre parameter read) | ! Special variable initialization | initialize_special ! (post parameter read) | ! Special variable finalization | finalize_special ! (deallocation, etc.) | ! | ! Special initial condition | init_special ! this is called last so may be used to modify | ! the mvar variables declared by this module | ! or optionally modify any of the other f array | ! variables. The latter, however, should be | ! avoided where ever possible. | ! | ! Special term in the mass (density) equation | special_calc_density ! Special term in the momentum (hydro) equation | special_calc_hydro ! Special term in the energy equation | special_calc_energy ! Special term in the induction (magnetic) | special_calc_magnetic ! equation | ! | ! Special equation | dspecial_dt ! NOT IMPLEMENTED FULLY YET - HOOKS NOT PLACED INTO THE PENCIL-CODE ! !** AUTOMATIC CPARAM.INC GENERATION **************************** ! Declare (for generation of special_dummies.inc) the number of f array ! variables and auxiliary variables added by this module ! ! CPARAM logical, parameter :: lspecial = .true. ! ! MVAR CONTRIBUTION 0 ! MAUX CONTRIBUTION 0 ! ! PENCILS PROVIDED uSH2; uRR2; uEL2; uRRm; bSH2; bRR2; bEL2; bRRm ! !** AUTOMATIC REFERENCE-LINK.TEX GENERATION ******************** ! Declare relevant citations from pencil-code/doc/citations/ref.bib for this module. ! The entries are taken from pencil-code/doc/citations/notes.tex ! !*************************************************************** ! ! HOW TO USE THIS FILE ! -------------------- ! ! Change the line above to ! lspecial = .true. ! to enable use of special hooks. ! ! The rest of this file may be used as a template for your own ! special module. Lines which are double commented are intended ! as examples of code. Simply fill out the prototypes for the ! features you want to use. ! ! Save the file with a meaningful name, eg. geo_kws.f90 and place ! it in the $PENCIL_HOME/src/special directory. This path has ! been created to allow users ot optionally check their contributions ! in to the Pencil-Code SVN repository. This may be useful if you ! are working on/using the additional physics with somebodyelse or ! may require some assistance from one of the main Pencil-Code team. ! ! To use your additional physics code edit the Makefile.local in ! the src directory under the run directory in which you wish to ! use your additional physics. Add a line with all the module ! selections to say something like: ! ! SPECIAL=special/geo_kws ! ! Where geo_kws it replaced by the filename of your new module ! upto and not including the .f90 ! module Special ! use Cdata use General, only: keep_compiler_quiet use Messages, only: svn_id, fatal_error ! implicit none ! include '../special.h' ! ! Declare index of new variables in f array (if any). ! logical :: luschur2_as_aux=.false., lbschur2_as_aux=.false. logical :: luschurm_as_aux=.false., lbschurm_as_aux=.false. logical :: luschurp2_as_aux=.false., lbschurp2_as_aux=.false. logical :: loee_as_aux=.false., leig_as_aux=.false. logical :: lovSH_as_aux=.false., lovRR_as_aux=.false. integer :: iuschur2, iuschur2_SH, iuschur2_RR, iuschur2_EL, iuschurm_RR integer :: ibschur2, ibschur2_SH, ibschur2_RR, ibschur2_EL, ibschurm_RR integer :: iuschurp2, iuschurp2_SH, iuschurp2_RR, iuschurp2_EL integer :: ibschurp2, ibschurp2_SH, ibschurp2_RR, ibschurp2_EL integer :: ioee, ieig, iovSH, iovRR ! logical :: luschur_as_aux=.false., lbschur_as_aux=.false. logical :: luse_complex_schur=.false. integer :: iuschur_SH, iuschur_RR, iuschur_EL integer :: ibschur_SH, ibschur_RR, ibschur_EL ! namelist /special_init_pars/ & luschur2_as_aux, lbschur2_as_aux, & luschurm_as_aux, lbschurm_as_aux, & luschurp2_as_aux, lbschurp2_as_aux, & luschur_as_aux, lbschur_as_aux, & loee_as_aux, leig_as_aux, & lovSH_as_aux, lovRR_as_aux ! logical :: luij_schur=.false., lbij_schur=.false., ldiagnos_always=.false. logical :: luschur_unprojected=.false., lbschur_unprojected=.false., lQ_schur_QT=.true. logical :: leigvec=.false. !PAR_DOC: compute eigenvectors ! namelist /special_run_pars/ & luij_schur, lbij_schur, ldiagnos_always, & luschur_unprojected, lbschur_unprojected, lQ_schur_QT, & luschur_as_aux, lbschur_as_aux, luse_complex_schur, leigvec, & loee_as_aux, leig_as_aux, & lovSH_as_aux, lovRR_as_aux ! ! Diagnostic variables (needs to be consistent with reset list below). ! integer :: idiag_uSH2=0 ! DIAG_DOC: $\left<{u^\mathrm{SH}}^2\right>$ integer :: idiag_uRR2=0 ! DIAG_DOC: $\left<{u^\mathrm{RR}}^2\right>$ integer :: idiag_uEL2=0 ! DIAG_DOC: $\left<{u^\mathrm{EL}}^2\right>$ integer :: idiag_uRRm=0 ! DIAG_DOC: $\left<2u^\mathrm{SH}u^\mathrm{RR}\right>$ integer :: idiag_bSH2=0 ! DIAG_DOC: $\left<{b^\mathrm{SH}}^2\right>$ integer :: idiag_bRR2=0 ! DIAG_DOC: $\left<{b^\mathrm{RR}}^2\right>$ integer :: idiag_bEL2=0 ! DIAG_DOC: $\left<{b^\mathrm{EL}}^2\right>$ integer :: idiag_bRRm=0 ! DIAG_DOC: $\left<2b^\mathrm{SH}b^\mathrm{RR}\right>$ ! contains !================== External SELECT function for DGEES ================== logical function selct(WR, WI) ! real, intent(in) :: WR, WI ! Return .TRUE. for eigenvalues to move to the TOP-LEFT block. ! We choose REAL eigenvalues => complex pair (if any) goes bottom-right. selct = (WI == 0.) ! end function selct !================== External SELECT function for ZGEES ================== ! Dummy SELECT function (not used because SORT = 'N') logical function select_dummy(w) ! complex, intent(in) :: w select_dummy = .false. ! end function select_dummy !**************************************************************************** subroutine register_special ! ! Set up indices for variables in special modules. ! ! 6-oct-03/tony: coded ! use FArrayManager use Sub, only: register_report_aux use SharedVariables, only: put_shared_variable ! if (lroot) call svn_id( & "$Id$") ! if (luschur2_as_aux) & call register_report_aux('uschur2', iuschur2, iuschur2_SH, iuschur2_RR, iuschur2_EL) ! if (lbschur2_as_aux) & call register_report_aux('bschur2', ibschur2, ibschur2_SH, ibschur2_RR, ibschur2_EL) ! if (luschurp2_as_aux) & call register_report_aux('uschurp2', iuschurp2, iuschurp2_SH, iuschurp2_RR, iuschurp2_EL) ! if (lbschurp2_as_aux) & call register_report_aux('bschurp2', ibschurp2, ibschurp2_SH, ibschurp2_RR, ibschurp2_EL) ! if (luschurm_as_aux) & call register_report_aux('uschurm_RR', iuschurm_RR) ! if (lbschurm_as_aux) & call register_report_aux('bschurm_RR', ibschurm_RR) ! if (luschur_as_aux) then call farray_register_auxiliary('uschur_SH', iuschur_SH, vector=9) call farray_register_auxiliary('uschur_RR', iuschur_RR, vector=9) call farray_register_auxiliary('uschur_EL', iuschur_EL, vector=9) endif ! if (lbschur_as_aux) then call farray_register_auxiliary('bschur_SH', ibschur_SH, vector=9) call farray_register_auxiliary('bschur_RR', ibschur_RR, vector=9) call farray_register_auxiliary('bschur_EL', ibschur_EL, vector=9) endif ! if (loee_as_aux) & call farray_register_auxiliary('oee', ioee, vector=3) ! if (leig_as_aux) & call farray_register_auxiliary('eig', ieig, vector=3) ! if (lovSH_as_aux) & call farray_register_auxiliary('ovSH', iovSH, vector=1) ! if (lovRR_as_aux) & call farray_register_auxiliary('ovRR', iovRR, vector=1) ! endsubroutine register_special !*********************************************************************** subroutine initialize_special(f) ! ! Called after reading parameters, but before the time loop. ! ! 06-oct-03/tony: coded ! use SharedVariables, only: get_shared_variable ! real, dimension (mx,my,mz,mfarray) :: f ! call keep_compiler_quiet(f) ! endsubroutine initialize_special !*********************************************************************** subroutine init_special(f) ! ! initialise special condition; called from start.f90 ! 06-oct-2003/tony: coded ! real, dimension (mx,my,mz,mfarray) :: f integer :: j ! intent(inout) :: f ! call keep_compiler_quiet(f) ! endsubroutine init_special !*********************************************************************** subroutine pencil_criteria_special ! if (loee_as_aux) lpenc_requested(i_oo)=.true. ! endsubroutine pencil_criteria_special !*********************************************************************** subroutine calc_diagnostic_pencils(f,p) ! ! Calculate Special pencils. ! Most basic pencils should come first, as others may depend on them. ! ! 24-nov-04/tony: coded ! 27-jan-26/axel: Since -r42600, redefined mixed term as the sum of mixed and original term ! 21-mar-26/TP: Moved here since do not seem to be needed on the rhs and thus do not need to be translated ! to the GPU ! use Sub, only: grad use General, only: eigvec3 ! real, dimension (mx,my,mz,mfarray) :: f type (pencil_case) :: p ! real, dimension (nx,3,3) :: SH, RR, EL real, dimension (3,3) :: matA, matV_SH, matV_RR, matV_EL, vv real, dimension (3) :: lambdas, vSH, vRR real, allocatable :: matV(:,:), matQ(:,:) complex, allocatable :: matV_cmplx(:,:) complex :: matQ_cmplx(3,3) real :: matA2=0., o2, oe, e2, vSH2, ovSH, vRR2, ovRR integer :: i,j,kk,ll, ij, l, nnn=3, ieigvec ! intent(inout) :: f intent(inout) :: p ! ! ! ! Possibility of applying triple decomposition of uij ! Do only when output onto command line. ! Redefine mixed term as the sum of mixed and original term. ! Complex Schur decomposition works only for %uSH2(l), p%uRR2(l), p%uEL2(l). ! if ((luij_schur .and. lout) .or. ldiagnos_always) then if (luse_complex_schur) then allocate(matV_cmplx(nnn,nnn)) else allocate(matV(nnn,nnn), matQ(nnn,nnn)) endif ! do l=1,nx matA=p%uij(l,:,:) matA2=sum(matA**2) if (matA2/=0.) then if (luse_complex_schur) then call schur_standardized_complex(matA, matV_cmplx, matQ_cmplx, p%uSH2(l), p%uRR2(l), p%uEL2(l)) else call schur_standardized(matA, matV, matQ) call schur_decompose(matV, matV_SH, matV_RR, matV_EL) p%uSH2(l)=sum(matV_SH**2) p%uRR2(l)=sum(matV_RR**2) p%uEL2(l)=sum(matV_EL**2) p%uRRm(l)=2.*sum(matV_SH*matV_RR)+p%uRR2(l) endif ! ! Compute o2 in several cases ! if (loee_as_aux .or. & lovSH_as_aux .or. & lovRR_as_aux) o2=p%oo(l,1)**2+p%oo(l,2)**2+p%oo(l,3)**2 ! ! Possibility of computing the three eigenvectors of the symmetric part of A ! and their alignment with the vorticity vector. ! if (leigvec) then call eigvec3(matA, vv, lambdas) if (loee_as_aux) then do ieigvec=1,3 e2=vv(1,ieigvec)**2+vv(2,ieigvec)**2+vv(3,ieigvec)**2 oe=p%oo(l,1)*vv(1,ieigvec)+p%oo(l,2)*vv(2,ieigvec)+p%oo(l,3)*vv(3,ieigvec) f(l1+l-1,m,n,ioee+ieigvec-1)=oe/(tini+sqrt(o2*e2)) enddo endif ! ! Output of eigenvalues. ! if (leig_as_aux) then do ieigvec=1,3 f(l1+l-1,m,n,ieig+ieigvec-1)=lambdas(ieigvec) enddo endif endif endif ! ! Possibility of uSH, uRR, and uEL matrices as auxiliary arrays ! Decomposition in the form GradU = Q (SH+RR+EL) Q^T, so ! GradU_ij = Q_ik (SH+RR+EL)_kl QT_lj = Q_ik Q_jl (SH+RR+EL)_kl. ! Setting lQ_schur_QT=.true. applies to GradU = Q^T (SH+RR+EL) Q. ! Otherwise, we'd have GradU = Q^T (SH+RR+EL) Q. ! if (luschur_as_aux .or. luschurp2_as_aux) then if (luschur_unprojected) then do i=1,3 do j=1,3 ij=3*(i-1)+(j-1) if (luse_complex_schur) then f(l1+l-1,m,n,iuschur_SH+ij)= abs(matV_cmplx(i,j)) f(l1+l-1,m,n,iuschur_RR+ij)= real(matV_cmplx(i,j)) f(l1+l-1,m,n,iuschur_EL+ij)=aimag(matV_cmplx(i,j)) else f(l1+l-1,m,n,iuschur_SH+ij)=matV_SH(i,j) f(l1+l-1,m,n,iuschur_RR+ij)=matV_RR(i,j) f(l1+l-1,m,n,iuschur_EL+ij)=matV_EL(i,j) endif enddo enddo else SH=0. RR=0. EL=0. do i=1,3 do j=1,3 do kk=1,3 do ll=1,3 if (lQ_schur_QT) then SH(l,i,j)=SH(l,i,j)+matQ(i,kk)*matQ(j,ll)*matV_SH(kk,ll) RR(l,i,j)=RR(l,i,j)+matQ(i,kk)*matQ(j,ll)*matV_RR(kk,ll) EL(l,i,j)=EL(l,i,j)+matQ(i,kk)*matQ(j,ll)*matV_EL(kk,ll) else SH(l,i,j)=SH(l,i,j)+matV_SH(kk,ll)*matQ(kk,i)*matQ(ll,j) RR(l,i,j)=RR(l,i,j)+matV_RR(kk,ll)*matQ(kk,i)*matQ(ll,j) EL(l,i,j)=EL(l,i,j)+matV_EL(kk,ll)*matQ(kk,i)*matQ(ll,j) endif enddo enddo ! ! Option of auxiliary out. ! if (luschur_as_aux) then ij=3*(i-1)+(j-1) f(l1+l-1,m,n,iuschur_SH+ij)=SH(l,i,j) f(l1+l-1,m,n,iuschur_RR+ij)=RR(l,i,j) f(l1+l-1,m,n,iuschur_EL+ij)=EL(l,i,j) endif enddo enddo ! ! Output norm of projected matrix. ! if (luschurp2_as_aux) then f(l1+l-1,m,n,iuschurp2_SH)=sum(SH(l,:,:)**2) f(l1+l-1,m,n,iuschurp2_RR)=sum(RR(l,:,:)**2) f(l1+l-1,m,n,iuschurp2_EL)=sum(EL(l,:,:)**2) endif ! ! Vectors related to the antisymmetric matrices SH and RR. ! if (lovSH_as_aux) then vSH(1)=-.5*(SH(l,2,3)-SH(l,3,2)) vSH(2)=-.5*(SH(l,3,1)-SH(l,1,3)) vSH(3)=-.5*(SH(l,1,2)-SH(l,2,1)) vSH2=vSH(1)**2+vSH(2)**2+vSH(3)**2 ovSH=p%oo(l,1)*vSH(1)+p%oo(l,2)*vSH(2)+p%oo(l,3)*vSH(3) f(l1+l-1,m,n,iovSH)=ovSH/(tini+sqrt(o2*vSH2)) endif ! if (lovSH_as_aux) then vRR(1)=-.5*(RR(l,2,3)-RR(l,3,2)) vRR(2)=-.5*(RR(l,3,1)-RR(l,1,3)) vRR(3)=-.5*(RR(l,1,2)-RR(l,2,1)) vRR2=vRR(1)**2+vRR(2)**2+vRR(3)**2 ovRR=p%oo(l,1)*vRR(1)+p%oo(l,2)*vRR(2)+p%oo(l,3)*vRR(3) f(l1+l-1,m,n,iovRR)=ovRR/(tini+sqrt(o2*vRR2)) endif ! endif endif enddo if (luse_complex_schur) then deallocate(matV_cmplx) else deallocate(matV, matQ) endif ! ! Possibility of uSH2, uRR2, and uEL2 as auxiliary arrays ! if (luschur2_as_aux) then f(l1:l2,m,n,iuschur2_SH)=p%uSH2 f(l1:l2,m,n,iuschur2_RR)=p%uRR2 f(l1:l2,m,n,iuschur2_EL)=p%uEL2 endif ! ! Mixed term, but note that there is no p%uRRm with luse_complex_schur=T ! if (luschurm_as_aux) then if (luse_complex_schur) then call fatal_error('calc_pencils_special','no p%uRRm with luse_complex_schur=T') else f(l1:l2,m,n,iuschurm_RR)=p%uRRm endif endif endif ! ! Possibility of applying triple decomposition of bij ! Do only when output onto command line. ! Redefine mixed term as the sum of mixed and original term. ! if ((lbij_schur .and. lout) .or. ldiagnos_always) then if (luse_complex_schur) then allocate(matV_cmplx(nnn,nnn)) ! else ! allocate(matV(nnn,nnn), matQ(nnn,nnn)) endif do l=1,nx matA=p%bij(l,:,:) matA2=sum(matA**2) if (matA2/=0.) then allocate(matV(nnn,nnn), matQ(nnn,nnn)) if (luse_complex_schur) then call schur_standardized_complex(matA, matV_cmplx, matQ_cmplx, p%bSH2(l), p%bRR2(l), p%bEL2(l)) else call schur_standardized(matA,matV,matQ) call schur_decompose(matV, matV_SH, matV_RR, matV_EL) p%bSH2(l)=sum(matV_SH**2) p%bRR2(l)=sum(matV_RR**2) p%bEL2(l)=sum(matV_EL**2) p%bRRm(l)=2.*sum(matV_SH*matV_RR)+p%bRR2(l) endif deallocate(matV, matQ) endif ! ! Possibility of bSH, bRR, and bEL matrices as auxiliary arrays ! if (lbschur_as_aux .or. lbschurp2_as_aux) then if (lbschur_unprojected) then do i=1,3 do j=1,3 ij=3*(i-1)+(j-1) if (luse_complex_schur) then f(l1+l-1,m,n,ibschur_SH+ij)= abs(matV_cmplx(i,j)) f(l1+l-1,m,n,ibschur_RR+ij)= real(matV_cmplx(i,j)) f(l1+l-1,m,n,ibschur_EL+ij)=aimag(matV_cmplx(i,j)) else f(l1+l-1,m,n,ibschur_SH+ij)=matV_SH(i,j) f(l1+l-1,m,n,ibschur_RR+ij)=matV_RR(i,j) f(l1+l-1,m,n,ibschur_EL+ij)=matV_EL(i,j) endif enddo enddo else SH=0. RR=0. EL=0. do i=1,3 do j=1,3 do kk=1,3 do ll=1,3 if (lQ_schur_QT) then SH(l,i,j)=SH(l,i,j)+matQ(i,kk)*matQ(j,ll)*matV_SH(kk,ll) RR(l,i,j)=RR(l,i,j)+matQ(i,kk)*matQ(j,ll)*matV_RR(kk,ll) EL(l,i,j)=EL(l,i,j)+matQ(i,kk)*matQ(j,ll)*matV_EL(kk,ll) else SH(l,i,j)=SH(l,i,j)+matV_SH(kk,ll)*matQ(kk,i)*matQ(ll,j) RR(l,i,j)=RR(l,i,j)+matV_RR(kk,ll)*matQ(kk,i)*matQ(ll,j) EL(l,i,j)=EL(l,i,j)+matV_EL(kk,ll)*matQ(kk,i)*matQ(ll,j) endif enddo enddo if (lbschur_as_aux) then ij=3*(i-1)+(j-1) f(l1+l-1,m,n,ibschur_SH+ij)=SH(l,i,j) f(l1+l-1,m,n,ibschur_RR+ij)=RR(l,i,j) f(l1+l-1,m,n,ibschur_EL+ij)=EL(l,i,j) endif enddo enddo ! ! output norm of projected matrix ! if (lbschurp2_as_aux) then f(l1+l-1,m,n,ibschurp2_SH)=sum(SH(l,:,:)**2) f(l1+l-1,m,n,ibschurp2_RR)=sum(RR(l,:,:)**2) f(l1+l-1,m,n,ibschurp2_EL)=sum(EL(l,:,:)**2) endif endif endif enddo if (luse_complex_schur) then deallocate(matV_cmplx) ! else ! deallocate(matV, matQ) endif ! ! Possibility of bSH2, bRR2, and bEL2 as auxiliary arrays ! Use f(1,1,1,ibschur2_EL) to encode the correct update time. ! if (lbschur2_as_aux) then f(l1:l2,m,n,ibschur2_SH)=p%bSH2 f(l1:l2,m,n,ibschur2_RR)=p%bRR2 f(l1:l2,m,n,ibschur2_EL)=p%bEL2 endif ! ! Mixed term, but note that there is no p%uRRm with luse_complex_schur=T ! if (lbschurm_as_aux) then if (luse_complex_schur) then call fatal_error('calc_pencils_special','no p%bRRm with luse_complex_schur=T') else f(l1:l2,m,n,ibschurm_RR)=p%bRRm endif endif endif ! endsubroutine !*********************************************************************** subroutine calc_pencils_special(f,p) ! ! Calculate Special pencils. ! Most basic pencils should come first, as others may depend on them. ! ! 24-nov-04/tony: coded ! 27-jan-26/axel: Since -r42600, redefined mixed term as the sum of mixed and original term ! use Sub, only: grad use General, only: eigvec3 ! real, dimension (mx,my,mz,mfarray) :: f type (pencil_case) :: p ! intent(inout) :: f intent(inout) :: p call calc_diagnostic_pencils(f,p) endsubroutine calc_pencils_special !*********************************************************************** subroutine schur_decompose(matV, matV_SH, matV_RR, matV_EL) ! real, dimension (3,3) :: matV_SH, matV_RR, matV_EL, matV_rest, matV_rest_trans real, allocatable :: matV(:,:) integer :: i, j, nnn=3 ! ! Shear: ! matV_SH(1,1)=0. matV_SH(1,2)=matV(1,2) matV_SH(1,3)=matV(1,3) ! matV_SH(2,1)=0. matV_SH(2,2)=0. matV_SH(2,3)=sign(1.,matV(2,3))*max(abs(matV(2,3))-abs(matV(3,2)),0.) ! matV_SH(3,1)=0. matV_SH(3,2)=sign(1.,matV(3,2))*max(abs(matV(3,2))-abs(matV(2,3)),0.) matV_SH(3,3)=0. ! ! Rest: ! matV_rest=matV-matV_SH ! ! Transposed matrix ! do j=1,3 do i=1,3 matV_rest_trans(i,j)=matV_rest(j,i) enddo enddo ! ! Antisymmetric and symmetric parts ! matV_RR=.5*(matV_rest-matV_rest_trans) matV_EL=.5*(matV_rest+matV_rest_trans) ! endsubroutine schur_decompose !*********************************************************************** subroutine print_mat(M) real , intent(in) :: M(:,:) integer :: r, c do r = 1, size(M,1) write(*,'(100(1X,F12.6))') ( M(r,c), c=1,size(M,2) ) enddo end subroutine print_mat !*********************************************************************** subroutine schur_standardized(A_input, T, VS) integer :: nnn, lda, ldvs, info, sdim, lwork integer :: i, j real :: tol real :: A_input(3,3) real , allocatable :: A(:,:), T(:,:), VS(:,:), WR(:), WI(:), WORK(:) logical , allocatable :: BWORK(:) real :: a11, a12, a21, a22 real :: rt1r, rt1i, rt2r, rt2i, cs, sn real :: R11, R12, R21, R22 real , allocatable :: TMP(:,:), TMP2(:,:) external dgees, dlanv2 ! LAPACK externals ! ---- Example matrix (3x3). Replace with your data. ---- nnn = 3 lda = nnn ldvs = nnn tol = 1.0d-12 allocate(A(nnn,nnn), WR(nnn), WI(nnn)) allocate(BWORK(nnn)) A = reshape(A_input, shape(A_input), order=(/2,1/) ) T = A ! ---- Workspace query for DGEES ---- lwork = -1 allocate(WORK(1)) call dgees('V','S', selct, nnn, T, lda, sdim, WR, WI, VS, ldvs, WORK, lwork, BWORK, info) if (info .ne. 0) then print *, 'DGEES(work query) failed, INFO=', info stop 1 end if lwork = int(WORK(1)) deallocate(WORK) allocate(WORK(lwork)) ! ---- Real Schur with sorting: real eigenvalues first (complex pair -> bottom-right) ---- T = A call dgees('V','S', selct, nnn, T, lda, sdim, WR, WI, VS, ldvs, WORK, lwork, BWORK, info) if (info .ne. 0) then print *, 'DGEES failed, INFO=', info stop 1 end if ! ---- Standardize the lower-right 2x2 block if present ---- if (nnn >= 2) then i = nnn-1 j = nnn if (abs(T(j,i)) > tol) then ! We have a 2x2 block at (i:i+1, i:i+1) a11 = T(i,i); a12 = T(i,j) a21 = T(j,i); a22 = T(j,j) ! DLANV2 computes rotation (cs,sn) that brings the 2x2 to standard real Schur form call dlanv2(a11, a12, a21, a22, rt1r, rt1i, rt2r, rt2i, cs, sn) ! Build R = [ cs sn; -sn cs ] R11 = cs; R12 = sn R21 = -sn; R22 = cs ! Apply similarity on T: T := (I ⊕ R^T) * T * (I ⊕ R) allocate(TMP(nnn,2), TMP2(2,nnn)) ! Right-multiply columns i:i+1 by R TMP(:,1) = T(:,i)*R11 + T(:,j)*R21 TMP(:,2) = T(:,i)*R12 + T(:,j)*R22 T(:,i) = TMP(:,1) T(:,j) = TMP(:,2) ! Left-multiply rows i:i+1 by R^T TMP2(1,:) = R11*T(i,:) + R12*T(j,:) TMP2(2,:) = R21*T(i,:) + R22*T(j,:) T(i,:) = TMP2(1,:) T(j,:) = TMP2(2,:) ! Update VS (Q): VS := VS * (I ⊕ R) TMP(:,1) = VS(:,i)*R11 + VS(:,j)*R21 TMP(:,2) = VS(:,i)*R12 + VS(:,j)*R22 VS(:,i) = TMP(:,1) VS(:,j) = TMP(:,2) deallocate(TMP, TMP2) ! Enforce upper-right >= 0 (beta >= 0) by a column/row sign flip if needed if (T(i,j) < 0.0d0) then ! Multiply by F = diag(1, -1): columns T(:,j) = -T(:,j) VS(:,j) = -VS(:,j) ! and rows T(j,:) = -T(j,:) end if ! Clean tiny roundoff if (abs(T(j,i)) < tol) T(j,i) = 0.0d0 if (abs(T(i,j) + T(j,i)) < tol) T(i,j) = -T(j,i) end if end if ! ---- Print results ---- !print *, 'T (standardized Schur form):' !call print_mat(T) !print *, 'Q (VS from DGEES):' !call print_mat(VS) endsubroutine schur_standardized !*********************************************************************** subroutine schur_standardized_complex(A_input, T, VS, SH2, RR2, EL2) integer :: nnn, lda, ldvs, info, sdim, lwork, i, j real :: A_input(3,3), rwork(3), SH2, RR2, EL2 real, allocatable :: A(:,:), WI(:), WORK(:) complex, allocatable :: T(:,:) complex :: W(3), VS(3,3) logical :: BWORK(1) ! SORT='N', so dummy size is sufficient external zgees ! LAPACK externals ! ---- Example matrix (3x3). Replace with your data. ---- nnn = 3 lda = nnn ldvs = nnn allocate(A(nnn,nnn), WI(nnn)) A = reshape(A_input, shape(A_input), order=(/2,1/) ) T = cmplx(A) lwork = 12 allocate(WORK(lwork)) ! ---- Real Schur with sorting: real eigenvalues first (complex pair -> bottom-right) ---- T = cmplx(A) call zgees('V','S', select_dummy, nnn, T, lda, sdim, W, VS, 3, WORK, lwork, rwork, BWORK, info) deallocate(WORK) ! SH2= abs(T(1,2))**2+ abs(T(1,3))**2+ abs(T(2,3))**2 RR2=aimag(T(1,1))**2+aimag(T(2,2))**2+aimag(T(3,3))**2 EL2= real(T(1,1))**2+ real(T(2,2))**2+ real(T(3,3))**2 ! if (info .ne. 0) then print *, 'ZGEES failed, INFO=', info stop 1 end if endsubroutine schur_standardized_complex !*********************************************************************** subroutine calc_diagnostics_special(f,p) use Diagnostics, only: sum_mn_name real, dimension (mx,my,mz,mfarray) :: f type (pencil_case) :: p ! intent(in) :: f,p if (ldiagnos .or. ldiagnos_always) then call sum_mn_name(p%uSH2,idiag_uSH2) call sum_mn_name(p%uRR2,idiag_uRR2) call sum_mn_name(p%uEL2,idiag_uEL2) call sum_mn_name(p%uRRm,idiag_uRRm) call sum_mn_name(p%bSH2,idiag_bSH2) call sum_mn_name(p%bRR2,idiag_bRR2) call sum_mn_name(p%bEL2,idiag_bEL2) call sum_mn_name(p%bRRm,idiag_bRRm) endif endsubroutine calc_diagnostics_special !*********************************************************************** subroutine dspecial_dt(f,df,p) ! ! The entire module could be renamed to Klein-Gordon or Scalar field equation. ! Calculate right hand side of ONE OR MORE extra coupled PDEs ! along the 'current' Pencil, i.e. f(l1:l2,m,n) where ! m,n are global variables looped over in equ.f90 ! ! Due to the multi-step Runge Kutta timestepping used one MUST always ! add to the present contents of the df array. NEVER reset it to zero. ! ! Several precalculated Pencils of information are passed for ! efficiency. ! ! 6-oct-03/tony: coded ! 2-nov-21/axel: first set of equations coded ! ! real, dimension (mx,my,mz,mfarray) :: f real, dimension (mx,my,mz,mvar) :: df type (pencil_case) :: p ! intent(in) :: f,p intent(inout) :: df ! ! Identify module and boundary conditions. ! if (headtt.or.ldebug) print*,'dspecial_dt: SOLVE dspecial_dt' call calc_diagnostics_special(f,p) ! endsubroutine dspecial_dt !*********************************************************************** subroutine dspecial_dt_ode ! endsubroutine dspecial_dt_ode !*********************************************************************** subroutine read_special_init_pars(iostat) ! use File_io, only: parallel_unit ! integer, intent(out) :: iostat ! read(parallel_unit, NML=special_init_pars, IOSTAT=iostat) ! endsubroutine read_special_init_pars !*********************************************************************** subroutine write_special_init_pars(unit) ! integer, intent(in) :: unit ! write(unit, NML=special_init_pars) ! endsubroutine write_special_init_pars !*********************************************************************** subroutine read_special_run_pars(iostat) ! use File_io, only: parallel_unit ! integer, intent(out) :: iostat ! read(parallel_unit, NML=special_run_pars, IOSTAT=iostat) ! endsubroutine read_special_run_pars !*********************************************************************** subroutine write_special_run_pars(unit) ! integer, intent(in) :: unit ! write(unit, NML=special_run_pars) ! endsubroutine write_special_run_pars !*********************************************************************** subroutine rprint_special(lreset,lwrite) ! ! Reads and registers print parameters relevant to special. ! ! 06-oct-03/tony: coded ! use Diagnostics, only: parse_name ! integer :: iname logical :: lreset,lwrite ! ! check for those quantities for which we want video slices ! if (lwrite_slices) then where(cnamev=='uschur2_SH' .or. cnamev=='uschur2_RR' .or. cnamev=='uschur2_EL' .or. & cnamev=='uschurm_RR' .or. & cnamev=='bschur2_SH' .or. cnamev=='bschur2_RR' .or. cnamev=='bschur2_EL' .or. & cnamev=='bschurm_RR' .or. & cnamev=='uschurp2_SH'.or. cnamev=='uschurp2_RR'.or. cnamev=='uschurp2_EL'.or. & cnamev=='bschurp2_SH'.or. cnamev=='bschurp2_RR'.or. cnamev=='bschurp2_EL' ) cformv='DEFINED' endif ! ! reset everything in case of reset ! (this needs to be consistent with what is defined above!) ! if (lreset) then idiag_uSH2=0; idiag_uRR2=0; idiag_uEL2=0; idiag_uRRm=0 idiag_bSH2=0; idiag_bRR2=0; idiag_bEL2=0; idiag_bRRm=0 endif ! do iname=1,nname call parse_name(iname,cname(iname),cform(iname),'uSH2',idiag_uSH2) call parse_name(iname,cname(iname),cform(iname),'uRR2',idiag_uRR2) call parse_name(iname,cname(iname),cform(iname),'uEL2',idiag_uEL2) call parse_name(iname,cname(iname),cform(iname),'uRRm',idiag_uRRm) call parse_name(iname,cname(iname),cform(iname),'bSH2',idiag_bSH2) call parse_name(iname,cname(iname),cform(iname),'bRR2',idiag_bRR2) call parse_name(iname,cname(iname),cform(iname),'bEL2',idiag_bEL2) call parse_name(iname,cname(iname),cform(iname),'bRRm',idiag_bRRm) enddo ! endsubroutine rprint_special !*********************************************************************** subroutine get_slices_special(f,slices) ! ! Write slices for animation of Special variables. ! ! 27-jan-26/axel: adapted from src/special/chiral_mhd.f90 ! use Slices_methods, only: assign_slices_scal ! real, dimension (mx,my,mz,mfarray) :: f type (slice_data) :: slices ! ! Loop over slices. ! select case (trim(slices%name)) case ('uschur2_SH'); call assign_slices_scal(slices,f,iuschur2_SH) case ('uschur2_RR'); call assign_slices_scal(slices,f,iuschur2_RR) case ('uschur2_EL'); call assign_slices_scal(slices,f,iuschur2_EL) case ('uschurm_RR'); call assign_slices_scal(slices,f,iuschurm_RR) case ('bschur2_SH'); call assign_slices_scal(slices,f,ibschur2_SH) case ('bschur2_RR'); call assign_slices_scal(slices,f,ibschur2_RR) case ('bschur2_EL'); call assign_slices_scal(slices,f,ibschur2_EL) case ('bschurm_RR'); call assign_slices_scal(slices,f,ibschurm_RR) case ('uschurp2_SH'); call assign_slices_scal(slices,f,iuschurp2_SH) case ('uschurp2_RR'); call assign_slices_scal(slices,f,iuschurp2_RR) case ('uschurp2_EL'); call assign_slices_scal(slices,f,iuschurp2_EL) case ('bschurp2_SH'); call assign_slices_scal(slices,f,ibschurp2_SH) case ('bschurp2_RR'); call assign_slices_scal(slices,f,ibschurp2_RR) case ('bschurp2_EL'); call assign_slices_scal(slices,f,ibschurp2_EL) endselect ! endsubroutine get_slices_special !*********************************************************************** subroutine special_after_boundary(f) ! ! Possibility to modify the f array after the boundaries are ! communicated. ! ! 16-apr-25/axel: coded ! ! use Mpicomm, only: mpireduce_sum, mpiallreduce_sum, mpibcast_real ! use Sub, only: dot2_mn, grad, curl, dot_mn ! real, dimension (mx,my,mz,mfarray), intent(in) :: f ! endsubroutine special_after_boundary !******************************************************************** !******************************************************************** !************ DO NOT DELETE THE FOLLOWING ************* !******************************************************************** !** This is an automatically generated include file that creates ** !** copies dummy routines from nospecial.f90 for any Special ** !** routines not implemented in this file ** !** ** include '../special_dummies.inc' !*********************************************************************** endmodule Special