Lambda_solve Subroutine

public subroutine Lambda_solve(spos_in, hmap_in, hmap_xv, X1_base_in, X2_base_in, LA_fbase_in, X1_in, X2_in, LA_s, phiPrime_s, chiPrime_s)

Uses

  • proc~~lambda_solve~~UsesGraph proc~lambda_solve Lambda_solve module~modgvec_base MODgvec_base proc~lambda_solve->module~modgvec_base module~modgvec_fbase MODgvec_fBase proc~lambda_solve->module~modgvec_fbase module~modgvec_globals MODgvec_Globals proc~lambda_solve->module~modgvec_globals module~modgvec_hmap MODgvec_hmap proc~lambda_solve->module~modgvec_hmap module~modgvec_base->module~modgvec_fbase module~modgvec_base->module~modgvec_globals module~modgvec_sbase MODgvec_sBase module~modgvec_base->module~modgvec_sbase module~modgvec_sgrid MODgvec_sGrid module~modgvec_base->module~modgvec_sgrid module~modgvec_fbase->module~modgvec_globals iso_fortran_env iso_fortran_env module~modgvec_globals->iso_fortran_env 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_c_hmap->module~modgvec_globals module~modgvec_hmap_axisnb->module~modgvec_fbase module~modgvec_hmap_axisnb->module~modgvec_globals module~modgvec_hmap_axisnb->module~modgvec_c_hmap module~modgvec_io_netcdf MODgvec_IO_NETCDF 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_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~modgvec_sgrid->module~modgvec_globals module~modgvec_io_netcdf->module~modgvec_globals 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

Solve for lambda on one given flux surface (spos_in), using weak form of J^s=0: d/dzeta(B_theta)-d/dtheta(B_zeta)=0

Note that the mapping defined by X1 and X2 must be fully initialized, since derivatives in s must be taken!

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in) :: spos_in

s position to evaluate lambda
class(c_hmap), intent(in) :: hmap_in
class(c_hmap_auxvar), intent(in) :: hmap_xv(X1_base_in%f%mn_IP)
class(t_base), intent(in) :: X1_base_in
class(t_base), intent(in) :: X2_base_in
type(t_fBase), intent(in) :: LA_fbase_in
real(kind=wp), intent(in) :: X1_in(1:X1_base_in%s%nBase,1:X1_base_in%f%modes)

U%X1 variable, is reshaped to 2D at input
real(kind=wp), intent(in) :: X2_in(1:X2_base_in%s%nBase,1:X2_base_in%f%modes)

U%X2 variable, is reshaped to 2D at input
real(kind=wp), intent(out) :: LA_s(1:LA_fbase_in%modes)

lambda at spos
real(kind=wp), intent(in) :: phiPrime_s

toroidal flux derivative phi' at the position s
real(kind=wp), intent(in) :: chiPrime_s

poloidal flux derivative chi' at the position s

Calls

proc~~lambda_solve~~CallsGraph proc~lambda_solve Lambda_solve proc~fbase_evaldof_ip_tens t_fBase%fBase_evalDOF_IP_tens proc~lambda_solve->proc~fbase_evaldof_ip_tens proc~hmap_eval_gij_aux c_hmap%hmap_eval_gij_aux proc~lambda_solve->proc~hmap_eval_gij_aux proc~hmap_eval_jh_aux c_hmap%hmap_eval_Jh_aux proc~lambda_solve->proc~hmap_eval_jh_aux proc~lambda_setup_and_solve Lambda_setup_and_solve proc~lambda_solve->proc~lambda_setup_and_solve proc~sbase_evaldof_s t_sBase%sBase_evalDOF_s proc~lambda_solve->proc~sbase_evaldof_s 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 eval_gij eval_gij proc~hmap_eval_gij_aux->eval_gij eval_Jh eval_Jh proc~hmap_eval_jh_aux->eval_Jh proc~fbase_projectiptodof_tens t_fBase%fBase_projectIPtoDOF_tens proc~lambda_setup_and_solve->proc~fbase_projectiptodof_tens proc~solve SOLVE proc~lambda_setup_and_solve->proc~solve proc~sbase_eval t_sBase%sBase_eval proc~sbase_evaldof_s->proc~sbase_eval proc~sbase_evaldof_base t_sBase%sBase_evalDOF_base proc~sbase_evaldof_s->proc~sbase_evaldof_base dgemv dgemv proc~fbase_evaldof_xn->dgemv proc~fbase_eval_xn t_fBase%fBase_eval_xn proc~fbase_evaldof_xn->proc~fbase_eval_xn proc~fbase_projectiptodof_tens->dgemm eval_basis eval_basis proc~sbase_eval->eval_basis eval_basis_and_n_derivs eval_basis_and_n_derivs proc~sbase_eval->eval_basis_and_n_derivs lagrangeinterpolationpolys lagrangeinterpolationpolys proc~sbase_eval->lagrangeinterpolationpolys proc~sgrid_find_elem t_sGrid%sGrid_find_elem proc~sbase_eval->proc~sgrid_find_elem dgetrf dgetrf proc~solve->dgetrf dgetrs dgetrs proc~solve->dgetrs

Called by

proc~~lambda_solve~~CalledByGraph proc~lambda_solve Lambda_solve proc~init_la_from_solution Init_LA_from_Solution proc~init_la_from_solution->proc~lambda_solve proc~initsolutionmhd3d t_functional_mhd3d%InitSolutionMHD3D proc~initsolutionmhd3d->proc~init_la_from_solution proc~initsolution InitSolution proc~initsolution->proc~initsolutionmhd3d proc~rungvec rungvec proc~rungvec->proc~initsolutionmhd3d proc~start_rungvec start_rungvec proc~start_rungvec->proc~rungvec program~gvec GVEC program~gvec->proc~rungvec

Source Code

SUBROUTINE Lambda_solve(spos_in,hmap_in,hmap_xv,X1_base_in,X2_base_in,LA_fbase_in,X1_in,X2_in,LA_s,phiPrime_s,chiPrime_s)
! MODULES
  USE MODgvec_Globals,       ONLY:n_warnings_occured
  USE MODgvec_base          ,ONLY: t_base
  USE MODgvec_fbase         ,ONLY: t_fbase
  USE MODgvec_hmap          ,ONLY: PP_T_HMAP,PP_T_HMAP_AUXVAR
  IMPLICIT NONE
!-----------------------------------------------------------------------------------------------------------------------------------
! INPUT VARIABLES
  CLASS(t_base),INTENT(IN)  :: X1_base_in,X2_base_in           !< base classes belong to solution X1_in,X2_in
  TYPE(t_fbase),INTENT(IN) :: LA_fbase_in                     !< base class belong to solution LA_s
#ifdef PP_WHICH_HMAP
  TYPE(PP_T_HMAP), INTENT(IN) :: hmap_in
  TYPE(PP_T_HMAP_AUXVAR), INTENT(IN) :: hmap_xv(X1_base_in%f%mn_IP)  !< auxiliary variables for hmap, must be pre-computed
#else
  CLASS(PP_T_HMAP), INTENT(IN) :: hmap_in
  CLASS(PP_T_HMAP_AUXVAR), INTENT(IN) :: hmap_xv(X1_base_in%f%mn_IP)  !< auxiliary variables for hmap, must be pre-computed
#endif
  REAL(wp)     , INTENT(IN) :: spos_in                  !! s position to evaluate lambda
  REAL(wp)     , INTENT(IN) :: X1_in(1:X1_base_in%s%nBase,1:X1_base_in%f%modes) !! U%X1 variable, is reshaped to 2D at input
  REAL(wp)     , INTENT(IN) :: X2_in(1:X2_base_in%s%nBase,1:X2_base_in%f%modes) !! U%X2 variable, is reshaped to 2D at input
  REAL(wp)     , INTENT(IN) :: phiPrime_s !! toroidal flux derivative phi' at the position s
  REAL(wp)     , INTENT(IN) :: chiPrime_s !! poloidal flux derivative chi' at the position s
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
REAL(wp)     , INTENT(  OUT) :: LA_s(1:LA_fbase_in%modes) !! lambda at spos
!-----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
  INTEGER                               :: iMode,i_mn,mn_IP
  REAL(wp)                              :: spos,minJ
  REAL(wp),DIMENSION(1:X1_base_in%f%modes) :: X1_s,X1_ds !! X1 solution at spos
  REAL(wp),DIMENSION(1:X2_base_in%f%modes) :: X2_s,X2_ds !! X1 solution at spos
  REAL(wp),DIMENSION(1:X1_base_in%f%mn_IP) :: X1_s_IP,dX1ds,dX1dthet,dX1dzeta, & !mn_IP should be same for all!
                                              X2_s_IP,dX2ds,dX2dthet,dX2dzeta, &
                                              detJ,gam_tt,gam_tz,gam_zz

!===================================================================================================================================
  __PERFON('lambda_solve')

  spos=MIN(1.0_wp-1.0e-12_wp,MAX(1.0e-04_wp,spos_in))
  mn_IP = X1_base_in%f%mn_IP
  IF(X2_base_in%f%mn_IP.NE.mn_IP) CALL abort(__STAMP__,&
                                             'X2 mn_IP /= X1 mn_IP')
  IF(LA_fbase_in%mn_IP .NE.mn_IP)  CALL abort(__STAMP__,&
                                             'LA mn_IP /= X1 mn_IP')


!$OMP PARALLEL DO        &
!$OMP   SCHEDULE(STATIC) &
!$OMP   DEFAULT(NONE)    &
!$OMP   PRIVATE(iMode)   &
!$OMP   SHARED(spos,X1_s,X1_ds,X1_base_in,X1_in)
  DO iMode=1,X1_base_in%f%modes
    X1_s( iMode)  = X1_base_in%s%evalDOF_s(spos,      0,X1_in(:,iMode))
    X1_ds(iMode)  = X1_base_in%s%evalDOF_s(spos,DERIV_S,X1_in(:,iMode))
  END DO
!$OMP END PARALLEL DO

!$OMP PARALLEL DO        &
!$OMP   SCHEDULE(STATIC) &
!$OMP   DEFAULT(NONE)    &
!$OMP   PRIVATE(iMode)   &
!$OMP   SHARED(spos,X2_s,X2_ds,X2_base_in,X2_in)
  DO iMode=1,X2_base_in%f%modes
    X2_s( iMode)  = X2_base_in%s%evalDOF_s(spos,      0,X2_in(:,iMode))
    X2_ds(iMode)  = X2_base_in%s%evalDOF_s(spos,DERIV_S,X2_in(:,iMode))
  END DO
!$OMP END PARALLEL DO

  X1_s_IP  = X1_base_in%f%evalDOF_IP(         0,X1_s )
  dX1ds    = X1_base_in%f%evalDOF_IP(         0,X1_ds)
  dX1dthet = X1_base_in%f%evalDOF_IP(DERIV_THET,X1_s )
  dX1dzeta = X1_base_in%f%evalDOF_IP(DERIV_ZETA,X1_s )

  X2_s_IP  = X2_base_in%f%evalDOF_IP(         0,X2_s )
  dX2ds    = X2_base_in%f%evalDOF_IP(         0,X2_ds)
  dX2dthet = X2_base_in%f%evalDOF_IP(DERIV_THET,X2_s )
  dX2dzeta = X2_base_in%f%evalDOF_IP(DERIV_ZETA,X2_s )



!$OMP PARALLEL DO        &
!$OMP   SCHEDULE(STATIC) &
!$OMP   DEFAULT(NONE)    &
!$OMP   PRIVATE(i_mn)  &
!$OMP   SHARED(mn_IP,X1_s_IP,X2_s_IP,dX1ds,dX2ds,dX1dthet,dX2dthet,detJ,hmap_in,hmap_xv)
  DO i_mn=1,mn_IP
    detJ(i_mn)= (dX1ds(i_mn)*dX2dthet(i_mn)-dX1dthet(i_mn)*dX2ds(i_mn)) &
               *hmap_in%eval_Jh_aux(X1_s_IP(i_mn),X2_s_IP(i_mn),hmap_xv(i_mn)) !J_p*J_h
  END DO !i_mn
!$OMP END PARALLEL DO

  minJ=MINVAL(detJ)
  IF(minJ.LT.1.0e-12) THEN
    n_warnings_occured=n_warnings_occured+1
    i_mn= MINLOC(detJ(:),1)
    WRITE(UNIT_stdOut,'(4X,A8,I8,4(A,E11.3))')'WARNING ',n_warnings_occured, &
                                                 ' : min(J)= ',MINVAL(detJ),' at s= ',spos, &
                                                                       ' theta= ',X1_base_in%f%x_IP(1,i_mn), &
                                                                        ' zeta= ',X1_base_in%f%x_IP(2,i_mn)
    i_mn= MAXLOC(detJ(:),1)
    WRITE(UNIT_stdOut,'(4X,16X,4(A,E11.3))')'     ...max(J)= ',MAXVAL(detJ),' at s= ',spos, &
                                                                       ' theta= ',X1_base_in%f%x_IP(1,i_mn), &
                                                                        ' zeta= ',X1_base_in%f%x_IP(2,i_mn)
!    CALL abort(__STAMP__, &
!        'Lambda_solve: Jacobian smaller that  1.0e-12!!!' )
  END IF

!$OMP PARALLEL DO        &
!$OMP   SCHEDULE(STATIC) &
!$OMP   DEFAULT(NONE)    &
!$OMP   PRIVATE(i_mn)  &
!$OMP   SHARED(mn_IP,gam_tt,gam_tz,gam_zz,detJ,hmap_in,hmap_xv,X1_s_IP,X2_s_IP,dX1dthet,dX2dthet,dX1dzeta,dX2dzeta)
  DO i_mn=1,mn_IP
    gam_tt(i_mn) = hmap_in%eval_gij_aux(dX1dthet(i_mn),dX2dthet(i_mn),0.0_wp, &
                                         X1_s_IP(i_mn), X2_s_IP(i_mn),        &
                                        dX1dthet(i_mn),dX2dthet(i_mn),0.0_wp, &
                                        hmap_xv(i_mn)) / detJ(i_mn)
    gam_tz(i_mn) = hmap_in%eval_gij_aux(dX1dthet(i_mn),dX2dthet(i_mn),0.0_wp, &
                                         X1_s_IP(i_mn), X2_s_IP(i_mn),        &
                                        dX1dzeta(i_mn),dX2dzeta(i_mn),1.0_wp, &
                                        hmap_xv(i_mn)) / detJ(i_mn)
    gam_zz(i_mn) = hmap_in%eval_gij_aux(dX1dzeta(i_mn),dX2dzeta(i_mn),1.0_wp, &
                                         X1_s_IP(i_mn), X2_s_IP(i_mn),        &
                                        dX1dzeta(i_mn),dX2dzeta(i_mn),1.0_wp, &
                                        hmap_xv(i_mn)) / detJ(i_mn)
  END DO !i_mn
!$OMP END PARALLEL DO

  CALL lambda_setup_and_solve(LA_fbase_in,phiPrime_s,ChiPrime_s,gam_tt,gam_tz,gam_zz,LA_s)

  __PERFOFF('lambda_solve')

END SUBROUTINE Lambda_solve