module ohm_electromagnetic_braginskii_m !! This module solves the electromagnetic and inertial part in Ohm's law !! !! \f[ !! \beta\frac{\partial}{\partial t}A_\parallel+\frac{\partial}{\partial t}\frac{j_\parallel}{n}=rhs, !! \f] !! with Ampere's law !! \f[ !! \nabla_\perp^2A_\parallel = -j_\parallel !! \f] use NETCDF use comm_handler_m, only : comm_handler_t use precision_grillix_m, only : GP, GP_NAN use error_handling_grillix_m, only: handle_error_netcdf use multistep_m, only : karniadakis_t use parallel_map_m, only : parallel_map_t use penalisation_m, only : penalisation_t use boundaries_braginskii_m, only : boundaries_braginskii_t ! From PARALLAX use perf_m, only : perf_start, perf_stop use screen_io_m, only : get_stdout use equilibrium_m, only : equilibrium_t use mesh_cart_m, only : mesh_cart_t use boundaries_perp_m, only : set_perpbnds, extrapolate_ghost_points use helmholtz_operator_m, only : helmholtz_single_inner use helmholtz_solver_m, only : helmholtz_solver_t use helmholtz_netcdfio_m, only : write_netcdf_helmholtz ! Parameters use params_brag_model_m, only : & rhos, beta, mass_ratio_ei, etapar0 use params_brag_switches_m, only : swb_ohm_physdiss use params_brag_boundaries_perp_m, only : & bnddescr_apar_core, bnddescr_apar_wall, bnddescr_apar_dome, bnddescr_apar_out use params_brag_boundaries_parpen_m, only : & pen_epsinv implicit none public :: ohm_advance public :: compute_jpar public :: apar_solve contains subroutine ohm_advance(comm_handler, & equi, mesh_stag, hsolver_stag, map, penalisation_stag, & boundaries, tstep_ohm, & nev_stag, nev_adv_stag, tev_adv_stag, & aparv, jparv, & dohm, & apar_adv, jpar_adv, & sinfo, res) !! Advances Ohm's law in time and solves for parallel elctromagnetic potential and parallel current type(comm_handler_t), intent(in) :: comm_handler !! Communicators class(equilibrium_t), intent(inout) :: equi !! Equilibrium type(mesh_cart_t), intent(in) :: mesh_stag !! Mesh (staggered) class(helmholtz_solver_t), intent(inout) :: hsolver_stag !! Elliptic (2D) solver on staggered mesh type(parallel_map_t), intent(in) :: map !! Parallel map type(penalisation_t), intent(in) :: penalisation_stag !! Penalisation (staggered) type(boundaries_braginskii_t), intent(inout) :: boundaries !! Boundary information for the BRAGINSKII model type(karniadakis_t), intent(inout) :: tstep_ohm !! Time-step integrator for Ohm's law real(GP), dimension(mesh_stag%get_n_points()), intent(in) :: nev_stag !! Values of density on staggered grid at time t real(GP), dimension(mesh_stag%get_n_points()), intent(in) :: nev_adv_stag !! Values of density on staggered grid at time t+1 real(GP), dimension(mesh_stag%get_n_points()), intent(in) :: tev_adv_stag !! Values of electron temperature on staggered grid at time t+1 real(GP), dimension(mesh_stag%get_n_points()), intent(in) :: aparv !! Values of parallel electromagnetic potential at time t real(GP), dimension(mesh_stag%get_n_points()), intent(in) :: jparv !! Values of parallel current at time t real(GP), dimension(mesh_stag%get_n_points()), intent(in) :: dohm !! Right hand side of Ohm's law real(GP), dimension(mesh_stag%get_n_points()), intent(inout) :: apar_adv !! Values of parallel electromagnetic potential at time t+1, on input values for initial guess real(GP), dimension(mesh_stag%get_n_points()), intent(out) :: jpar_adv !! Values of parallel current at time t+1 integer, intent(out) :: sinfo !! Info from solver real(GP), intent(out) :: res !! Residual of solution real(GP), dimension(mesh_stag%get_n_points()) :: psipar, psipar_adv !! Generalised parallel electromagnetic potential at t and t+1 integer :: ki, l call perf_start('../../ohm_advance') ! Compute generalised electromagnetic potential call perf_start('../../../compute_psipar') !$omp parallel default(none) private(l) & !$omp shared(mesh_stag, beta, mass_ratio_ei, nev_stag, aparv, jparv, psipar) !$omp do do l = 1, mesh_stag%get_n_points() psipar(l) = beta * aparv(l) + mass_ratio_ei * jparv(l) / nev_stag(l) ! note: values on ghost points not actually needed in the following, here kept for convenience and simpler interface enddo !$omp end do !$omp end parallel ! Advance in time call tstep_ohm%advance(dohm, psipar, psipar_adv) call perf_stop('../../../compute_psipar') ! Solve for apar^t+1 call perf_start('../../../apar_solve') call apar_solve(comm_handler, equi, mesh_stag, hsolver_stag, & boundaries, & psipar_adv, apar_adv, sinfo, res, nev_adv_stag, & tev_adv_stag, penalisation_stag, tstep_ohm) call perf_stop('../../../apar_solve') ! Compute jpar^t+1 call perf_start('../../../compute_jpar') call compute_jpar(equi, mesh_stag, boundaries, apar_adv, jpar_adv) call perf_stop('../../../compute_jpar') call tstep_ohm%shift_storage((/psipar, dohm/)) call perf_stop('../../ohm_advance') end subroutine subroutine compute_jpar(equi, mesh_stag, boundaries, aparv, jparv) !! Computes parallel current from parallel electromagnetic potential !! \f[ !! j_\parallel = -\nabla_\perp^2 A_\parallel !! \f] class(equilibrium_t), intent(inout) :: equi !! Equilibrium type(mesh_cart_t), intent(in) :: mesh_stag !! Mesh within poloidal plane type(boundaries_braginskii_t), intent(inout) :: boundaries real(GP), dimension(mesh_stag%get_n_points()), intent(in) :: aparv !! Values of parallel electromagnetic potential real(GP), dimension(mesh_stag%get_n_points()), intent(out) :: jparv !! Values of parallel current integer :: ki, kb, kg, l real(GP) :: phi_stag, x, y, fac, xi real(GP), dimension(mesh_stag%get_n_points()) :: co ! Set up coefficients fac = rhos**2 phi_stag = mesh_stag%get_phi() !$omp parallel default(none) & !$omp private(ki, kb, kg, l, x, y, xi) & !$omp shared(equi, mesh_stag, phi_stag, fac, boundaries, & !$omp co, jparv, aparv) ! setup coefficients !$omp do do ki = 1, mesh_stag%get_n_points_inner() l = mesh_stag%inner_indices(ki) x = mesh_stag%get_x(l) y = mesh_stag%get_y(l) co(l) = equi%jacobian(x, y, phi_stag) enddo !$omp end do !$omp do do kb = 1, mesh_stag%get_n_points_boundary() l = mesh_stag%boundary_indices(kb) x = mesh_stag%get_x(l) y = mesh_stag%get_y(l) co(l) = equi%jacobian(x, y, phi_stag) enddo !$omp end do !$omp do do kg = 1, mesh_stag%get_n_points_ghost() l = mesh_stag%ghost_indices(kg) x = mesh_stag%get_x(l) y = mesh_stag%get_y(l) co(l) = equi%jacobian(x, y, phi_stag) enddo !$omp end do ! compute jpar on inner grid !$omp do do ki = 1, mesh_stag%get_n_points_inner() l = mesh_stag%inner_indices(ki) x = mesh_stag%get_x(l) y = mesh_stag%get_y(l) xi = fac / equi%jacobian(x, y, phi_stag) jparv(l) = helmholtz_single_inner(mesh_stag, aparv, l, co, xi) enddo !$omp end do ! Set boundaries call set_perpbnds(mesh_stag, boundaries%jpar%types, jparv, boundaries%jpar%vals) ! Extrapolate ghost points call extrapolate_ghost_points(mesh_stag, jparv) !$omp end parallel end subroutine subroutine apar_solve(comm_handler, equi, mesh_stag, hsolver_stag, boundaries, & psiparv, aparv, sinfo, res, co_inert, co_te, penalisation_stag, tstep_ohm) !! Solves for parallel electromagnetic potential !! if co_inert presents, \f[ !! \beta A_\parallel -\frac{mu}{n}\nabla_\perp^2A_\parallel= \psi_\parallel !! \f] !! otherwise, \f[ !! -\nabla_\perp^2A_\parallel= \psi_\parallel !! \f] type(comm_handler_t), intent(in) :: comm_handler !! Communicators class(equilibrium_t), intent(in) :: equi !! Equilibrium type(mesh_cart_t), intent(in) :: mesh_stag !! Mesh within poloidal plane class(helmholtz_solver_t), intent(inout) :: hsolver_stag !! Elliptic (2D) solver on staggered mesh type(boundaries_braginskii_t), intent(inout) :: boundaries !! Boundary information for the BRAGINSKII model real(GP), dimension(mesh_stag%get_n_points()), intent(in) :: psiparv !! Values for generalised electromagnetic potential real(GP), dimension(mesh_stag%get_n_points()), intent(inout) :: aparv !! Values for parallel electromagnetic potential, on input initial guess integer, intent(out) :: sinfo !! Info from solver real(GP), intent(out) :: res !! Residual of solution real(GP), dimension(mesh_stag%get_n_points()), intent(in), optional :: co_inert !! Coefficient in inertial term (= Values of density at time t+1 on staggered grid) real(GP), dimension(mesh_stag%get_n_points()), intent(in), optional :: co_te !! Coefficient in implicit resistivity (= Values of electron temperature at time t+1 on staggered grid) type(penalisation_t), intent(in), optional :: penalisation_stag !! Penalisation (staggered) type(karniadakis_t), intent(inout), optional :: tstep_ohm !! Time-step integrator for Ohm's law integer :: ki, kb,kg, l, nf90_stat, nf90_id real(GP) :: phi_stag, x, y, fac, pen_fac1, pen_fac2 real(GP), dimension(mesh_stag%get_n_points()) :: co, rhs real(GP), dimension(mesh_stag%get_n_points_inner()) :: lambda, xi character(len=3) :: plane_c fac = rhos**2 phi_stag = mesh_stag%get_phi() ! Setup equation system (coefficients and right hand side) call perf_start('../../../../ohm_setup') !$omp parallel default(none) & !$omp private(ki, kb, kg, l, x, y, pen_fac1, pen_fac2) & !$omp shared(equi, mesh_stag, phi_stag, boundaries, & !$omp beta, mass_ratio_ei, etapar0, fac, swb_ohm_physdiss, & !$omp lambda, xi, co, rhs, psiparv, & !$omp co_inert, co_te, penalisation_stag, pen_epsinv, tstep_ohm) if (present(co_inert).and.present(co_te).and.present(penalisation_stag).and.present(tstep_ohm)) then !$omp do do ki = 1, mesh_stag%get_n_points_inner() l = mesh_stag%inner_indices(ki) x = mesh_stag%get_x(l) y = mesh_stag%get_y(l) pen_fac1 = 1.0_GP + penalisation_stag%get_charfun_val(ki) * pen_epsinv * tstep_ohm%get_dtau_bdf() pen_fac2 = 1.0_GP - penalisation_stag%get_charfun_val(ki) co(l) = equi%jacobian(x, y, phi_stag) lambda(ki) = pen_fac1 * beta xi(ki) = pen_fac1 * fac * mass_ratio_ei / (co_inert(l) * equi%jacobian(x, y, phi_stag)) & + pen_fac2 * fac * tstep_ohm%get_dtau_bdf() * swb_ohm_physdiss * etapar0 & / (sqrt(co_te(l)**3) * equi%jacobian(x, y, phi_stag)) rhs(l) = psiparv(l) enddo !$omp end do else !$omp do do ki = 1, mesh_stag%get_n_points_inner() l = mesh_stag%inner_indices(ki) x = mesh_stag%get_x(l) y = mesh_stag%get_y(l) co(l) = equi%jacobian(x, y, phi_stag) lambda(ki) = 0.0_GP xi(ki) = fac / equi%jacobian(x, y, phi_stag) rhs(l) = psiparv(l) enddo !$omp end do endif !$omp do do kb = 1, mesh_stag%get_n_points_boundary() l = mesh_stag%boundary_indices(kb) x = mesh_stag%get_x(l) y = mesh_stag%get_y(l) co(l) = equi%jacobian(x, y, phi_stag) rhs(l) = boundaries%apar%vals(kb) enddo !$omp end do !$omp do do kg = 1, mesh_stag%get_n_points_ghost() l = mesh_stag%ghost_indices(kg) x = mesh_stag%get_x(l) y = mesh_stag%get_y(l) co(l) = equi%jacobian(x, y, phi_stag) rhs(l) = 0.0_GP enddo !$omp end do !$omp end parallel call perf_stop('../../../../ohm_setup') call perf_start('../../../../ohm_init') call hsolver_stag%update(co, lambda, xi, & bnddescr_apar_core, & bnddescr_apar_wall, & bnddescr_apar_dome, & bnddescr_apar_out) call perf_stop('../../../../ohm_init') call perf_start('../../../../ohm_solve') call hsolver_stag%solve(rhs, aparv, res, sinfo) call perf_stop('../../../../ohm_solve') if (sinfo < 0) then ! Write out state if solver failed write(plane_c,'(I3.3)')comm_handler%get_plane() ! Overwrites existing file nf90_stat = nf90_create('apar_solve_failstate_plane'//plane_c//'.nc', & NF90_NETCDF4+NF90_CLOBBER, nf90_id) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) call write_netcdf_helmholtz(nf90_id, mesh_stag, bnddescr_apar_core, & bnddescr_apar_wall, bnddescr_apar_dome, & bnddescr_apar_out, & co, lambda, xi, & rhs, & guess=aparv, & sol=aparv, & hcsr_write_on=.true.) nf90_stat = nf90_close(nf90_id) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) endif ! Extrapolate ghost points call perf_start('../../../../ohm_ghosts') !$omp parallel default(none) shared(mesh_stag, aparv) call extrapolate_ghost_points(mesh_stag, aparv) !$omp end parallel call perf_stop('../../../../ohm_ghosts') end subroutine end module