EvalAux Subroutine

public subroutine EvalAux(dofs_in, JacCheck)

Uses

  • proc~~evalaux~~UsesGraph proc~evalaux EvalAux module~modgvec_globals MODgvec_Globals proc~evalaux->module~modgvec_globals module~modgvec_mhd3d_vars MODgvec_MHD3D_Vars proc~evalaux->module~modgvec_mhd3d_vars module~modgvec_mpi MODgvec_MPI proc~evalaux->module~modgvec_mpi module~modgvec_sol_var_mhd3d MODgvec_sol_var_MHD3D proc~evalaux->module~modgvec_sol_var_mhd3d iso_fortran_env iso_fortran_env module~modgvec_globals->iso_fortran_env module~modgvec_mhd3d_vars->module~modgvec_globals module~modgvec_mhd3d_vars->module~modgvec_sol_var_mhd3d module~modgvec_base MODgvec_base module~modgvec_mhd3d_vars->module~modgvec_base module~modgvec_boundaryfromfile MODgvec_boundaryFromFile module~modgvec_mhd3d_vars->module~modgvec_boundaryfromfile module~modgvec_hmap MODgvec_hmap module~modgvec_mhd3d_vars->module~modgvec_hmap module~modgvec_rprofile_base MODgvec_rProfile_base module~modgvec_mhd3d_vars->module~modgvec_rprofile_base module~modgvec_sgrid MODgvec_sGrid module~modgvec_mhd3d_vars->module~modgvec_sgrid module~modgvec_sol_var_mhd3d->module~modgvec_globals module~modgvec_c_sol_var MODgvec_c_sol_var module~modgvec_sol_var_mhd3d->module~modgvec_c_sol_var module~modgvec_base->module~modgvec_globals module~modgvec_base->module~modgvec_sgrid module~modgvec_fbase MODgvec_fBase module~modgvec_base->module~modgvec_fbase module~modgvec_sbase MODgvec_sBase module~modgvec_base->module~modgvec_sbase module~modgvec_boundaryfromfile->module~modgvec_globals module~modgvec_io_netcdf MODgvec_IO_NETCDF module~modgvec_boundaryfromfile->module~modgvec_io_netcdf module~modgvec_c_sol_var->module~modgvec_globals module~modgvec_c_hmap MODgvec_c_hmap module~modgvec_hmap->module~modgvec_c_hmap module~modgvec_hmap_axisnb MODgvec_hmap_axisNB module~modgvec_hmap->module~modgvec_hmap_axisnb module~modgvec_hmap_cyl MODgvec_hmap_cyl module~modgvec_hmap->module~modgvec_hmap_cyl module~modgvec_hmap_frenet MODgvec_hmap_frenet module~modgvec_hmap->module~modgvec_hmap_frenet module~modgvec_hmap_knot MODgvec_hmap_knot module~modgvec_hmap->module~modgvec_hmap_knot module~modgvec_hmap_rz MODgvec_hmap_RZ module~modgvec_hmap->module~modgvec_hmap_rz module~modgvec_rprofile_base->module~modgvec_globals module~modgvec_sgrid->module~modgvec_globals module~modgvec_c_hmap->module~modgvec_globals module~modgvec_fbase->module~modgvec_globals module~modgvec_hmap_axisnb->module~modgvec_globals module~modgvec_hmap_axisnb->module~modgvec_c_hmap module~modgvec_hmap_axisnb->module~modgvec_fbase module~modgvec_hmap_axisnb->module~modgvec_io_netcdf module~modgvec_hmap_cyl->module~modgvec_globals module~modgvec_hmap_cyl->module~modgvec_c_hmap module~modgvec_hmap_frenet->module~modgvec_globals module~modgvec_hmap_frenet->module~modgvec_c_hmap module~modgvec_hmap_knot->module~modgvec_globals module~modgvec_hmap_knot->module~modgvec_c_hmap module~modgvec_hmap_rz->module~modgvec_globals module~modgvec_hmap_rz->module~modgvec_c_hmap module~modgvec_io_netcdf->module~modgvec_globals module~modgvec_sbase->module~modgvec_globals module~modgvec_sbase->module~modgvec_sgrid module~sll_m_bsplines sll_m_bsplines module~modgvec_sbase->module~sll_m_bsplines module~sll_m_spline_interpolator_1d sll_m_spline_interpolator_1d module~modgvec_sbase->module~sll_m_spline_interpolator_1d module~sll_m_spline_matrix sll_m_spline_matrix module~modgvec_sbase->module~sll_m_spline_matrix module~sll_m_assert sll_m_assert module~sll_m_bsplines->module~sll_m_assert module~sll_m_bsplines_base sll_m_bsplines_base module~sll_m_bsplines->module~sll_m_bsplines_base module~sll_m_bsplines_non_uniform sll_m_bsplines_non_uniform module~sll_m_bsplines->module~sll_m_bsplines_non_uniform module~sll_m_bsplines_uniform sll_m_bsplines_uniform module~sll_m_bsplines->module~sll_m_bsplines_uniform module~sll_m_errors sll_m_errors module~sll_m_bsplines->module~sll_m_errors module~sll_m_working_precision sll_m_working_precision module~sll_m_bsplines->module~sll_m_working_precision module~sll_m_spline_interpolator_1d->module~sll_m_spline_matrix module~sll_m_spline_interpolator_1d->module~sll_m_assert module~sll_m_boundary_condition_descriptors sll_m_boundary_condition_descriptors module~sll_m_spline_interpolator_1d->module~sll_m_boundary_condition_descriptors module~sll_m_spline_interpolator_1d->module~sll_m_bsplines_base module~sll_m_spline_interpolator_1d->module~sll_m_errors module~sll_m_spline_1d sll_m_spline_1d module~sll_m_spline_interpolator_1d->module~sll_m_spline_1d module~sll_m_spline_interpolator_1d->module~sll_m_working_precision module~sll_m_spline_matrix->module~sll_m_errors module~sll_m_spline_matrix_banded sll_m_spline_matrix_banded module~sll_m_spline_matrix->module~sll_m_spline_matrix_banded module~sll_m_spline_matrix_base sll_m_spline_matrix_base module~sll_m_spline_matrix->module~sll_m_spline_matrix_base module~sll_m_spline_matrix_dense sll_m_spline_matrix_dense module~sll_m_spline_matrix->module~sll_m_spline_matrix_dense module~sll_m_spline_matrix->module~sll_m_working_precision module~sll_m_boundary_condition_descriptors->module~sll_m_working_precision module~sll_m_bsplines_base->module~sll_m_assert module~sll_m_bsplines_base->module~sll_m_working_precision module~sll_m_bsplines_non_uniform->module~sll_m_assert module~sll_m_bsplines_non_uniform->module~sll_m_bsplines_base module~sll_m_bsplines_non_uniform->module~sll_m_working_precision module~sll_m_bsplines_uniform->module~sll_m_assert module~sll_m_bsplines_uniform->module~sll_m_bsplines_base module~sll_m_bsplines_uniform->module~sll_m_errors module~sll_m_bsplines_uniform->module~sll_m_working_precision module~sll_m_errors->iso_fortran_env module~sll_m_spline_1d->module~sll_m_assert module~sll_m_spline_1d->module~sll_m_bsplines_base module~sll_m_spline_1d->module~sll_m_working_precision module~sll_m_spline_matrix_banded->iso_fortran_env module~sll_m_spline_matrix_banded->module~sll_m_assert module~sll_m_spline_matrix_banded->module~sll_m_errors module~sll_m_spline_matrix_banded->module~sll_m_spline_matrix_base module~sll_m_spline_matrix_banded->module~sll_m_working_precision module~sll_m_spline_matrix_base->module~sll_m_working_precision module~sll_m_spline_matrix_dense->iso_fortran_env module~sll_m_spline_matrix_dense->module~sll_m_assert module~sll_m_spline_matrix_dense->module~sll_m_errors module~sll_m_spline_matrix_dense->module~sll_m_spline_matrix_base module~sll_m_spline_matrix_dense->module~sll_m_working_precision

Evaluate auxiliary variables at input state, writes onto module variables!!! This includes the check of the Jacobian, which is crucial for the validity of the computations. only a positive Jacobian is accepted.

Arguments

Type IntentOptional Attributes Name
class(t_sol_var_MHD3D), intent(in) :: dofs_in

input solution
integer, intent(inout) :: JacCheck

if 1 on input: abort if detJ<1e-12max|detJ|. if 2 on input, no abort, unchanged if min(detJ)>=1e-12max|detJ|, else return -3 negative Jacobian everywhere (left-handed coordinates), -2 relative Jacobian too small, -1 Jacobian with sign change.

Calls

proc~~evalaux~~CallsGraph proc~evalaux EvalAux eval_all eval_all proc~evalaux->eval_all interface~par_allreduce par_AllReduce proc~evalaux->interface~par_allreduce proc~base_evaldof t_base%base_evalDOF proc~evalaux->proc~base_evaldof proc~base_evaldof_all t_base%base_evalDOF_all proc~evalaux->proc~base_evaldof_all proc~par_allreduce_array1d par_AllReduce_array1D interface~par_allreduce->proc~par_allreduce_array1d proc~par_allreduce_array2d par_AllReduce_array2D interface~par_allreduce->proc~par_allreduce_array2d proc~par_allreduce_scalar par_AllReduce_scalar interface~par_allreduce->proc~par_allreduce_scalar proc~par_allreduce_scalar_int par_AllReduce_scalar_int interface~par_allreduce->proc~par_allreduce_scalar_int proc~fbase_evaldof_ip_tens t_fBase%fBase_evalDOF_IP_tens proc~base_evaldof->proc~fbase_evaldof_ip_tens proc~base_evaldof_all->proc~fbase_evaldof_ip_tens dgemm dgemm proc~fbase_evaldof_ip_tens->dgemm proc~fbase_evaldof_xn t_fBase%fBase_evalDOF_xn proc~fbase_evaldof_ip_tens->proc~fbase_evaldof_xn dgemv dgemv proc~fbase_evaldof_xn->dgemv proc~fbase_eval_xn t_fBase%fBase_eval_xn proc~fbase_evaldof_xn->proc~fbase_eval_xn

Called by

proc~~evalaux~~CalledByGraph proc~evalaux EvalAux proc~evalenergy EvalEnergy proc~evalenergy->proc~evalaux proc~evalforce EvalForce proc~evalforce->proc~evalaux proc~evaltotals EvalTotals proc~evaltotals->proc~evalaux proc~minimizemhd3d_resetdescent t_minimizer_mhd3d%MinimizeMHD3d_ResetDescent proc~minimizemhd3d_resetdescent->proc~evalaux proc~minimizemhd3d_resetdescent->proc~evalenergy proc~minimizemhd3d_resetdescent->proc~evalforce proc~initsolutionmhd3d t_functional_mhd3d%InitSolutionMHD3D proc~initsolutionmhd3d->proc~evalenergy proc~initsolutionmhd3d->proc~evalforce interface~writestate WriteState proc~initsolutionmhd3d->interface~writestate proc~minimizemhd3d_descent t_minimizer_mhd3d%MinimizeMHD3D_descent proc~minimizemhd3d_descent->proc~evalenergy proc~minimizemhd3d_descent->proc~evalforce proc~minimizemhd3d_descent->proc~minimizemhd3d_resetdescent proc~minimizemhd3d_descent->interface~writestate proc~writestatetoascii WriteStateToASCII proc~writestatetoascii->proc~evaltotals program~gvec_post GVEC_POST program~gvec_post->proc~evalenergy program~gvec_post->proc~evalforce interface~writestate->proc~writestatetoascii proc~initsolution~2 InitSolution proc~initsolution~2->proc~initsolutionmhd3d proc~minimizemhd3d t_functional_mhd3d%MinimizeMHD3D proc~minimizemhd3d->proc~minimizemhd3d_descent proc~rungvec rungvec proc~rungvec->proc~initsolutionmhd3d proc~rungvec->proc~minimizemhd3d proc~minimize minimize proc~minimize->proc~minimizemhd3d proc~start_rungvec start_rungvec proc~start_rungvec->proc~rungvec program~gvec GVEC program~gvec->proc~rungvec

Source Code

SUBROUTINE EvalAux(dofs_in,JacCheck)
! MODULES
  USE MODgvec_MPI             , ONLY: par_AllReduce
  USE MODgvec_Globals         , ONLY: n_warnings_occured,myRank
  USE MODgvec_MHD3D_vars      , ONLY: X1_base,X2_base,LA_base,hmap,hmap_auxvar
  USE MODgvec_sol_var_MHD3D   , ONLY: t_sol_var_MHD3D
  IMPLICIT NONE
!-----------------------------------------------------------------------------------------------------------------------------------
! INPUT VARIABLES
  CLASS(t_sol_var_MHD3D), INTENT(IN   ) :: dofs_in  !! input solution
  INTEGER               , INTENT(INOUT) :: JacCheck !! if 1 on input: abort if detJ<1e-12*max|detJ|.
                                                    !! if 2 on input, no abort, unchanged if min(detJ)>=1e-12*max|detJ|, else return
                                                    !! -3 negative Jacobian everywhere (left-handed coordinates),
                                                    !! -2 relative Jacobian too small,
                                                    !! -1 Jacobian with sign change.
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
!-----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
  INTEGER   :: iGP,i_mn,IP_GP(2),JacError
  REAL(wp)  :: min_detJ,max_detJ,min_abs_detJ,max_abs_detJ
!===================================================================================================================================
  IF(JacCheck.LE.-1) THEN
      CALL abort(__STAMP__, &
          'You already called EvalAux, with a bad Jacobian, detJ < 1e-12*max|detJ|!!!', IntInfo=JacCheck )
  END IF

  __PERFON('EvalAux')

  __PERFON('EvalDOF_1')
  !2D data: interpolation points x gauss-points
  !CALL X1_base%evalDOF((/0,0/)         ,Uin%X1,X1_IP_GP  )
  !CALL X1_base%evalDOF((/0,DERIV_THET/),Uin%X1,dX1_dthet )
  !CALL X1_base%evalDOF((/0,DERIV_ZETA/),Uin%X1,dX1_dzeta)
  CALL X1_base%evalDOF_all(dofs_in%X1, y_IP_GP=X1_IP_GP, &
                                   dy_dthet_IP_GP=dX1_dthet, &
                                   dy_dzeta_IP_GP=dX1_dzeta )
  CALL X1_base%evalDOF((/DERIV_S,0/)   , dofs_in%X1, dX1_ds    )

  !CALL X2_base%evalDOF((/0,0/)         ,Uin%X2,X2_IP_GP  )
  !CALL X2_base%evalDOF((/0,DERIV_THET/),Uin%X2,dX2_dthet )
  !CALL X2_base%evalDOF((/0,DERIV_ZETA/),Uin%X2,dX2_dzeta )
  CALL X2_base%evalDOF_all(dofs_in%X2, y_IP_GP=X2_IP_GP, &
                                  dy_dthet_IP_GP=dX2_dthet, &
                                   dy_dzeta_IP_GP=dX2_dzeta )
  CALL X2_base%evalDOF((/DERIV_S,0/)   , dofs_in%X2, dX2_ds    )

  __PERFOFF('EvalDOF_1')

  __PERFON('eval_hmap')
  CALL hmap%eval_all((/X1_base%f%mn_nyq(1),X1_base%f%mn_nyq(2),nGP_end-nGP_str+1/),2,hmap_auxvar, &
                     X1_IP_GP,X2_IP_GP,dX1_dthet,dX2_dthet,dX1_dzeta,dX2_dzeta, &
                     J_h,   g_tt, g_tz,g_zz, &
                     Jh_dq1,gtt_dq1,gtz_dq1,gzz_dq1, &
                     Jh_dq2,gtt_dq2,gtz_dq2,gzz_dq2, &
                     g_t1,g_t2,g_z1,g_z2,Gh11,Gh22)
  __PERFOFF('eval_hmap')

  __PERFON('loop_1')
  min_detJ =HUGE(1.0_wp)
  max_detJ =-HUGE(1.0_wp)
!$OMP PARALLEL DO        &
!$OMP   SCHEDULE(STATIC) DEFAULT(NONE)    &
!$OMP   PRIVATE(iGP,i_mn)  &
!$OMP   REDUCTION(min:min_detJ) &
!$OMP   REDUCTION(max:max_detJ) &
!$OMP   SHARED(nGP_str,nGP_end,mn_IP,J_p,J_h,detJ,dX1_ds,dX2_dthet,dX2_ds,dX1_dthet)
  DO iGP=nGP_str,nGP_end
    DO i_mn=1,mn_IP
      J_p(  i_mn,iGP) = ( dX1_ds(i_mn,iGP)*dX2_dthet(i_mn,iGP) &
                         -dX2_ds(i_mn,iGP)*dX1_dthet(i_mn,iGP) )

      detJ(i_mn,iGP) = J_p(i_mn,iGP)*J_h(i_mn,iGP)
      min_detJ  = MIN(min_detJ,detJ(i_mn,iGP))
      max_detJ  = MAX(max_detJ,detJ(i_mn,iGP))
    END DO !i_mn
  END DO !iGP
!$OMP END PARALLEL DO
  __PERFOFF('loop_1')
  CALL par_AllReduce(min_detJ,'MIN')
  CALL par_AllReduce(max_detJ,'MAX')
  !max|detJ| = max(|min(detJ)|,|max(detJ)|)
  max_abs_detJ = MAX(ABS(min_detJ),ABS(max_detJ))

  !check relative Jacobian J/max|J|
  JacError=0
  IF(min_detJ .LT.1.0e-12_wp*max_abs_detJ) THEN
    !detJ < 1e-12*max|detJ|
    IF(max_detJ .LT. -1.0e-12_wp*max_abs_detJ) THEN
      JacError=-3 ! all detJ is negative, left-handed coordinates found (max(detJ) < -1e-12*max|detJ|)
    ELSEIF(MIN(ABS(min_detJ),ABS(max_detJ)) .LE. 1.0e-12_wp*max_abs_detJ) THEN
      JacError=-2 ! relative Jacobian too small ( |min(detJ)|<= 1e-12*max|detJ| or |max(detJ)|<= 1e-12*max|detJ| )
    ELSE
      JacError=-1 ! negative Jacobian found with sign change  (min(detJ) < -1e-12*max|detJ| and max(detJ) > 1e-12*max|detJ| )
    END IF
  END IF
  IF(JacError.NE.0) THEN
    SELECT CASE(JacCheck)
    CASE(1)
      n_warnings_occured=n_warnings_occured+1
      IP_GP= MINLOC(detJ(:,nGP_str:nGP_end))
      WRITE(UNIT_stdOut,'(4X,A8,I8,4(A,E11.3))')'WARNING ',n_warnings_occured, &
           &                                       ' : min(J)= ',MINVAL(detJ(:,nGP_str:nGP_end)),' at s= ',s_GP(IP_GP(2)), &
           &                                                             ' theta= ',X1_base%f%x_IP(1,IP_GP(1)), &
           &                                                              ' zeta= ',X1_base%f%x_IP(2,IP_GP(1))
      IP_GP= MAXLOC(detJ(:,:))
      WRITE(UNIT_stdOut,'(4X,16X,4(A,E11.3))')'     ...max(J)= ',MAXVAL(detJ(:,nGP_str:nGP_end)),' at s= ',s_GP(IP_GP(2)), &
           &                                                             ' theta= ',X1_base%f%x_IP(1,IP_GP(1)), &
           &                                                              ' zeta= ',X1_base%f%x_IP(2,IP_GP(1))

      SELECT CASE(JacError)
      CASE(-3)
        CALL abort(__STAMP__, &
          'EvalAux: negative Jacobian found everywhere (max(detJ) < -1e-12*max|detJ|). Left-handed coordinates cannot be used!!!', IntInfo=myRank )
      CASE(-2)
        CALL abort(__STAMP__, &
          'EvalAux: relative Jacobian too small (|min(detJ)| <= 1e-12*max|detJ| or |max(detJ)|<= 1e-12*max|detJ|). !!!', IntInfo=myRank )
      CASE(-1)
        CALL abort(__STAMP__, &
          'EvalAux: negative Jacobian found with sign change (min(detJ) < -1e-12*max|detJ| and max(detJ) > 1e-12*max|detJ| )!!!', IntInfo=myRank )
      END SELECT

    CASE(2) !quiet check, give back
      JacCheck=JacError
    END SELECT
  ELSE
    JacCheck=1 !set to default for safety (abort if detJ<0)
  END IF

  IF(JacCheck.LE.-1) THEN
    __PERFOFF('EvalAux')
    RETURN
  END IF

  __PERFON('EvalDOF_2')
  !2D data: interpolation points x gauss-points
  !CALL LA_base%evalDOF((/0,DERIV_THET/),Uin%LA,dLA_dthet)
  !CALL LA_base%evalDOF((/0,DERIV_ZETA/),Uin%LA,dLA_dzeta)
  CALL LA_base%evalDOF_all(dofs_in%LA, dy_dthet_IP_GP=dLA_dthet, &
                                   dy_dzeta_IP_GP=dLA_dzeta )
  __PERFOFF('EvalDOF_2')


  __PERFON('loop_2')
!$OMP PARALLEL DO        &
!$OMP   SCHEDULE(STATIC) DEFAULT(NONE)    &
!$OMP   PRIVATE(iGP,i_mn)  &
!$OMP   SHARED(nGP_str,nGP_end,mn_IP,b_thet,b_zeta,sJ_p,sJ_h,sdetJ,sJ_bcov_thet,sJ_bcov_zeta,bbcov_sJ, &
!$OMP          J_p,J_h,chiPrime_GP,phiPrime_GP,dLA_dzeta,dLA_dthet,g_tt,g_tz,g_zz)
  DO iGP=nGP_str,nGP_end
    DO i_mn=1,mn_IP
      b_thet(i_mn,iGP) = (chiPrime_GP(iGP)- phiPrime_GP(iGP)*dLA_dzeta(i_mn,iGP))    !b_theta
      b_zeta(i_mn,iGP) = phiPrime_GP(iGP)*(1.0_wp          + dLA_dthet(i_mn,iGP))    !b_zeta
      sJ_p(i_mn,iGP)         = 1.0_wp/J_p(i_mn,iGP)
      sJ_h(i_mn,iGP)         = 1.0_wp/J_h(i_mn,iGP)
      sdetJ(i_mn,iGP)        =  sJ_p(i_mn,iGP)*  sJ_h(i_mn,iGP)
      sJ_bcov_thet(i_mn,iGP) = (g_tt(i_mn,iGP)*b_thet(i_mn,iGP) + g_tz(i_mn,iGP)*b_zeta(i_mn,iGP))*sdetJ(i_mn,iGP)
      sJ_bcov_zeta(i_mn,iGP) = (g_tz(i_mn,iGP)*b_thet(i_mn,iGP) + g_zz(i_mn,iGP)*b_zeta(i_mn,iGP))*sdetJ(i_mn,iGP)
      bbcov_sJ(    i_mn,iGP) =  b_thet( i_mn,iGP)*sJ_bcov_thet(i_mn,iGP) &
                               +b_zeta( i_mn,iGP)*sJ_bcov_zeta(i_mn,iGP)
    END DO !i_mn
  END DO !iGP
!$OMP END PARALLEL DO
  __PERFOFF('loop_2')

  __PERFOFF('EvalAux')

END SUBROUTINE EvalAux