Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

evaluation of moments improved #20

Merged
merged 4 commits into from
Nov 30, 2023
Merged

evaluation of moments improved #20

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
3d29516
moment evaluation improved
jonatanschatzlmayr Apr 21, 2023
12b6dad
started working on perturbations from mephit
jonatanschatzlmayr May 4, 2023
45325b7
fixed bug in circular mesh
jonatanschatzlmayr Nov 27, 2023
e91dcad
fixed time issue with moment integral
jonatanschatzlmayr Nov 28, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 5 additions & 0 deletions EXAMPLES/example_1/gorilla.inp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -38,4 +38,9 @@ boole_grid_for_find_tetra = .false. ,
a_factor = 4 ,
b_factor = 0 ,
c_factor = 0 ,
boole_save_electric = .false. ,
boole_strong_electric_field = .false. ,
filename_electric_drift = 'electric_drift.dat' ,
filename_electric_field = 'electric_field.dat' ,
boole_pert_from_mephit = .false. ,
/
5 changes: 5 additions & 0 deletions EXAMPLES/example_2/gorilla.inp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -38,4 +38,9 @@ boole_grid_for_find_tetra = .false. ,
a_factor = 4 ,
b_factor = 0 ,
c_factor = 0 ,
boole_save_electric = .false. ,
boole_strong_electric_field = .false. ,
filename_electric_drift = 'electric_drift.dat' ,
filename_electric_field = 'electric_field.dat' ,
boole_pert_from_mephit = .false. ,
/
5 changes: 5 additions & 0 deletions EXAMPLES/example_3/gorilla.inp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -38,4 +38,9 @@ a_factor = 4 ,
b_factor = 0 ,
c_factor = 0 ,
i_time_tracing_option = 1 ,
boole_save_electric = .false. ,
boole_strong_electric_field = .false. ,
filename_electric_drift = 'electric_drift.dat' ,
filename_electric_field = 'electric_field.dat' ,
boole_pert_from_mephit = .false. ,
/
5 changes: 5 additions & 0 deletions EXAMPLES/example_4/gorilla.inp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -38,4 +38,9 @@ boole_grid_for_find_tetra = .false. ,
a_factor = 4 ,
b_factor = 0 ,
c_factor = 0 ,
boole_save_electric = .false. ,
boole_strong_electric_field = .false. ,
filename_electric_drift = 'electric_drift.dat' ,
filename_electric_field = 'electric_field.dat' ,
boole_pert_from_mephit = .false. ,
/
1 change: 1 addition & 0 deletions EXAMPLES/example_5/gorilla.inp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -38,4 +38,5 @@
boole_strong_electric_field = .false. ,
filename_electric_drift = 'electric_drift.dat' ,
filename_electric_field = 'electric_field.dat' ,
boole_pert_from_mephit = .false. ,
/
5 changes: 5 additions & 0 deletions EXAMPLES/example_6/gorilla.inp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -38,4 +38,9 @@ a_factor = 4 ,
b_factor = 0 ,
c_factor = 0 ,
i_time_tracing_option = 1 ,
boole_save_electric = .false. ,
boole_strong_electric_field = .false. ,
filename_electric_drift = 'electric_drift.dat' ,
filename_electric_field = 'electric_field.dat' ,
boole_pert_from_mephit = .false. ,
/
5 changes: 5 additions & 0 deletions EXAMPLES/example_7/gorilla.inp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,4 +32,9 @@ non_axi_noise_eps_A = 0.0001 ,
poly_order = 2 ,
rel_err_ode45 = 1e-08 ,
i_time_tracing_option = 1 ,
boole_save_electric = .false. ,
boole_strong_electric_field = .false. ,
filename_electric_drift = 'electric_drift.dat' ,
filename_electric_field = 'electric_field.dat' ,
boole_pert_from_mephit = .false. ,
/
5 changes: 5 additions & 0 deletions INPUT/gorilla.inp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -177,6 +177,11 @@
filename_electric_drift = 'electric_drift.dat',
!
!------------------------------------------------------------------------------------------------------------!
!
!perturbations to the fields can optionally be delivered by Mephit
boole_pert_from_mephit = .false.
!
!------------------------------------------------------------------------------------------------------------!
!
/

Expand Down
8 changes: 8 additions & 0 deletions SRC/chamb_m.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,3 +1,9 @@
! module chamb_sub

! implicit none

! contains

!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!
subroutine chamb_can(y,phi,ierr)
Expand All @@ -22,3 +28,5 @@ subroutine chamb_can(y,phi,ierr)
endif
!
end subroutine chamb_can

!end module chamb_sub
10 changes: 5 additions & 5 deletions SRC/circular_mesh.f90
View file
Open in desktop
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The check was useless in its original position:

(a) the checked arrays were not allocated at this point, therefore size() would give any kind of output back
(b) it never threw an error nonetheless, as the .not. was not properly placed, so that the total output would actually only be true, if verts (what mostly the case due to (a)) did not align with n_tetra AND simultaneously, neighbour/neighbour_faces DID match (by some unlikely chance, as they were also not allocated before).

The position/check idea is properly an artifact of an earlies version of the code.

Original file line number Diff line number Diff line change
Expand Up @@ -243,11 +243,6 @@ subroutine calc_mesh(verts_per_ring, n_slices, points, n_tetras, &
n_verts = verts_per_slice * n_slices

n_tetras = calc_n_tetras(verts_per_ring, n_slices, repeat_center_point)
if (.not. size(verts, dim=2) == n_tetras .and. size(neighbours, dim=2) == n_tetras &
.and. size(neighbour_faces, dim=2) == n_tetras) then
print *, "Invalid size of output arrays for function calc_mesh"
stop
end if

tetras_per_slice = n_tetras / n_slices

Expand All @@ -266,6 +261,11 @@ subroutine calc_mesh(verts_per_ring, n_slices, points, n_tetras, &
perbou_phi(4, n_tetras), perbou_theta(4, n_tetras))
allocate(top_facing_prisms(verts_per_slice - 1))

if (.not. (size(verts, dim=2) == n_tetras .and. size(neighbours, dim=2) == n_tetras &
.and. size(neighbour_faces, dim=2) == n_tetras)) then
print *, "Invalid size of output arrays for function calc_mesh"
stop
end if
!
! 7------6
! /| /|
Expand Down
11 changes: 5 additions & 6 deletions SRC/circular_mesh_SOLEDGE3X_EIRENE.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -129,12 +129,6 @@ subroutine calc_mesh_SOLEDGE3X_EIRENE(n_slices, points_rphiz, verts_per_slice, n
n_verts = verts_per_slice * n_slices
!
n_tetras = n_triangles*3*n_slices
!
if (.not. size(verts, dim=2) == n_tetras .and. size(neighbours, dim=2) == n_tetras &
.and. size(neighbour_faces, dim=2) == n_tetras) then
print *, "Invalid size of output arrays for function calc_mesh_SOLEDGE3X_EIRENE"
stop
end if
!
tetras_per_slice = n_tetras / n_slices
!
Expand All @@ -146,6 +140,11 @@ subroutine calc_mesh_SOLEDGE3X_EIRENE(n_slices, points_rphiz, verts_per_slice, n
allocate(verts(4, n_tetras), neighbours(4, n_tetras), neighbour_faces(4, n_tetras), &
perbou_phi(4, n_tetras))
!
if (.not.(size(verts, dim=2) == n_tetras .and. size(neighbours, dim=2) == n_tetras &
.and. size(neighbour_faces, dim=2) == n_tetras)) then
print *, "Invalid size of output arrays for function calc_mesh_SOLEDGE3X_EIRENE"
stop
end if
!
! 7------6
! /| /|
Expand Down
6 changes: 3 additions & 3 deletions SRC/find_tetra_mod.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@ module find_tetra_mod
!
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!
subroutine grid_for_find_tetra
subroutine grid_for_find_tetra !compare with chapter 4.1 of master thesis of Jonatan Schatzlmayr
!
use constants, only: pi, eps
use tetra_grid_mod, only: verts_rphiz, verts_sthetaphi, tetra_grid, ntetr
Expand Down Expand Up @@ -83,12 +83,10 @@ subroutine grid_for_find_tetra
nc = 200*c_factor
boole_axi_symmetry = .true.
endif
!boole_axi_symmetry = .false.
!
delta_c = (cmax - cmin)/nc
!
if (b_factor.eq.0) b_factor = maxval((/nint(2*pi/(grid_size(2)*n_field_periods*sqrt((delta_a**2+delta_c**2)/2))),1/))
b_factor = b_factor*n_field_periods
!
nb = grid_size(2)*b_factor
delta_b = 2*pi/nb
Expand Down Expand Up @@ -442,6 +440,8 @@ subroutine find_tetra(x,vpar,vperp,ind_tetr_out,iface,sign_t_step_in)
end select
!
!print*, count(entry_counter.eq.0),size(entry_counter)
!
!search tetrahedra one by one
do i = 1,ntetr_searched
!
if (use_grid.eqv..true.) then
Expand Down
6 changes: 5 additions & 1 deletion SRC/gorilla_settings_mod.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -86,6 +86,9 @@ module gorilla_settings_mod
!Boolean for precalculation of rectangular grid to improve find_tetra (sensible for n_particles >> 1)
logical, public, protected :: boole_grid_for_find_tetra
integer, public :: a_factor, b_factor, c_factor
!
!Boolean for importing field perturbations from Mephit
logical, public, protected :: boole_pert_from_mephit
!
!Namelist for Gorilla input
NAMELIST /GORILLANML/ eps_Phi, coord_system, ispecies, ipusher, &
Expand All @@ -97,7 +100,8 @@ module gorilla_settings_mod
& boole_adaptive_time_steps, desired_delta_energy, max_n_intermediate_steps, boole_grid_for_find_tetra, &
& a_factor, b_factor, c_factor, &
& i_time_tracing_option, &
& boole_strong_electric_field, boole_save_electric, filename_electric_field, filename_electric_drift
& boole_strong_electric_field, boole_save_electric, filename_electric_field, filename_electric_drift, &
& boole_pert_from_mephit
!
public :: load_gorilla_inp,set_eps_Phi, optional_quantities_type
!
Expand Down
8 changes: 2 additions & 6 deletions SRC/orbit_timestep_gorilla.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -147,7 +147,7 @@ end subroutine orbit_timestep_gorilla
!
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!
subroutine initialize_gorilla(i_option,ipert_in,bmod_multiplier)!,boole_grid_for_find_tetra)
subroutine initialize_gorilla(i_option,ipert_in,bmod_multiplier)
!
use constants, only: echarge,ame,amp,clight
use gorilla_settings_mod, only: eps_Phi,coord_system,ispecies,boole_grid_for_find_tetra
Expand All @@ -161,7 +161,6 @@ subroutine initialize_gorilla(i_option,ipert_in,bmod_multiplier)!,boole_grid_for
implicit none
!
integer :: iper,ipert
!logical, intent(inout), optional :: boole_grid_for_find_tetra
integer,intent(in),optional :: i_option,ipert_in
logical :: boole_grid,boole_physics,boole_bmod_multiplier
double precision,intent(in),optional :: bmod_multiplier
Expand Down Expand Up @@ -263,10 +262,7 @@ subroutine initialize_gorilla(i_option,ipert_in,bmod_multiplier)!,boole_grid_for
!
call make_precomp_poly()
!
!if (present(boole_grid_for_find_tetra)) then
!if (grid_kind.eq.1) boole_grid_for_find_tetra = .false. !rectangular grid
if (boole_grid_for_find_tetra) call grid_for_find_tetra
!endif
if (boole_grid_for_find_tetra) call grid_for_find_tetra
!
endif !boole_physics
!
Expand Down
133 changes: 36 additions & 97 deletions SRC/pusher_tetra_poly.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2162,6 +2162,7 @@ subroutine calc_optional_quantities(poly_order,z0,tau,optional_quantities)
use gorilla_settings_mod, only: boole_time_Hamiltonian, boole_gyrophase, boole_vpar_int, boole_vpar2_int, &
& optional_quantities_type
use tetra_physics_mod, only: hamiltonian_time,cm_over_e,tetra_physics
use tetra_grid_mod, only: verts_rphiz, tetra_grid
!
implicit none
!
Expand All @@ -2178,121 +2179,59 @@ subroutine calc_optional_quantities(poly_order,z0,tau,optional_quantities)
double precision, dimension(:), allocatable :: vpar_coef,res_poly_coef
double precision, dimension(4,poly_order+1) :: dummy_coef_mat !exists because of analytic_coeff_without_precomp syntax
logical, dimension(4) :: boole_faces_off = .false. !avoid root coefficients computation
double precision, dimension(:), allocatable :: dtdtau_coef, omega_coef
integer :: i
!
!recalculate polynomial coefficients (tensors) if more than one integration step was performed
if (number_of_integration_steps .gt. 1) call set_integration_coef_manually(poly_order,z0)
!
allocate(x_coef(3,poly_order+1))
allocate(vpar_coef(poly_order+1))
!
call z_series_coef(poly_order,z0,x_coef,vpar_coef)
!
if(boole_time_hamiltonian.or.boole_vpar_int.or.boole_vpar2_int) then
allocate(x_coef(3,poly_order+1))
allocate(vpar_coef(poly_order+1))
allocate(dtdtau_coef(poly_order+1))
allocate(x_vpar_coef(3,poly_order+1))
!
allocate(x_vpar_coef(3,poly_order+1))
call z_series_coef(poly_order,z0,x_coef,vpar_coef)
x_vpar_coef = poly_multiplication(vpar_coef,x_coef)
!
call poly_multiplication_coef(x_coef(1,:),vpar_coef(:),x_vpar_coef(1,:))
call poly_multiplication_coef(x_coef(2,:),vpar_coef(:),x_vpar_coef(2,:))
call poly_multiplication_coef(x_coef(3,:),vpar_coef(:),x_vpar_coef(3,:))
endif
do i = 1,poly_order+1
dtdtau_coef(i) = sum(hamiltonian_time(ind_tetr)%vec_mismatch_der * x_coef(:,i)) + &
& cm_over_e * hamiltonian_time(ind_tetr)%h1_in_curlh * vpar_coef(i) + &
& cm_over_e * sum(hamiltonian_time(ind_tetr)%vec_parcurr_der * x_vpar_coef(:,i))
enddo
dtdtau_coef(1) = dtdtau_coef(1) + hamiltonian_time(ind_tetr)%h1_in_curlA
!
!Optional computation of Hamiltonian time
if(boole_time_hamiltonian) then
delta_t_hamiltonian = hamiltonian_time(ind_tetr)%h1_in_curlA * tau + &
& cm_over_e * hamiltonian_time(ind_tetr)%h1_in_curlh * moment_integration(poly_order,tau,vpar_coef)+&
& sum( hamiltonian_time(ind_tetr)%vec_mismatch_der * moment_integration(poly_order,tau,x_coef) ) + &
& cm_over_e * sum(hamiltonian_time(ind_tetr)%vec_parcurr_der * &
& moment_integration(poly_order,tau,x_vpar_coef))
!
!Multiply delta_t_hamiltonian with appropriate sign (We require that tau remains positive inside the algorithm)
delta_t_hamiltonian = delta_t_hamiltonian * dble(sign_rhs)
!
! call calc_t_hamiltonian(poly_order,z0,tau,delta_t_hamiltonian)
!
optional_quantities%t_hamiltonian = optional_quantities%t_hamiltonian + delta_t_hamiltonian
!
!Optional computation of gyrophase
if(boole_gyrophase) then
optional_quantities%gyrophase = optional_quantities%gyrophase - ( &
!
!Zeroth term
& 1.d0/cm_over_e * tetra_physics(ind_tetr)%bmod1 * delta_t_hamiltonian + &
!
!First term
& 1.d0/cm_over_e * sum( tetra_physics(ind_tetr)%gb * moment_integration(poly_order,tau,x_coef) ) *&
& hamiltonian_time(ind_tetr)%h1_in_curlA * dble(sign_rhs) + &
!
!Second term
& sum(moment_integration(poly_order,tau,x_vpar_coef) * tetra_physics(ind_tetr)%gb ) * &
& hamiltonian_time(ind_tetr)%h1_in_curlh * dble(sign_rhs)+ &
!
!Third term
& 1.d0/cm_over_e * sum( matmul( moment_integration(poly_order,tau,poly_multiplication(x_coef,x_coef)) , &
& hamiltonian_time(ind_tetr)%vec_mismatch_der ) * tetra_physics(ind_tetr)%gb) * dble(sign_rhs) + &
!
!Fourth term
& sum( matmul( moment_integration(poly_order,tau,poly_multiplication(x_coef,x_vpar_coef)), &
& hamiltonian_time(ind_tetr)%vec_parcurr_der ) * tetra_physics(ind_tetr)%gb) * dble(sign_rhs) &
& )
optional_quantities%t_hamiltonian = optional_quantities%t_hamiltonian + moment_integration(poly_order,tau,dtdtau_coef)&
&*dble(sign_rhs)
endif
!

endif ! gyrophase
!Optional computation of gyrophase
if(boole_gyrophase) then
allocate(omega_coef(poly_order+1))
do i = 1,poly_order+1
omega_coef(i) = 1.d0/cm_over_e * sum(tetra_physics(ind_tetr)%gb * x_coef(:,i))
enddo
omega_coef(1) = omega_coef(1) + 1.d0/cm_over_e * tetra_physics(ind_tetr)%bmod1
!
endif ! time_Hamiltonian
optional_quantities%gyrophase = optional_quantities%gyrophase - dble(sign_rhs) * &
& moment_integration(poly_order,tau,poly_multiplication(dtdtau_coef,omega_coef))
deallocate(omega_coef)
endif
!
!Optional computation of integral of v_par
if (boole_vpar_int) then
! optional_quantities%vpar_int = optional_quantities%vpar_int + moment_integration(poly_order,tau,vpar_coef)* &
! & dt_dtau_const
optional_quantities%vpar_int = optional_quantities%vpar_int + &
!
!First term
& hamiltonian_time(ind_tetr)%h1_in_curlA * moment_integration(poly_order,tau,vpar_coef) + &
!
!Second term
& cm_over_e * hamiltonian_time(ind_tetr)%h1_in_curlh * &
& moment_integration(poly_order,tau,poly_multiplication(vpar_coef,vpar_coef))+ &
!
!Third term
& sum( hamiltonian_time(ind_tetr)%vec_mismatch_der * &
& moment_integration(poly_order,tau,x_vpar_coef) ) + &
!
!Fourth term
& cm_over_e * sum(hamiltonian_time(ind_tetr)%vec_parcurr_der * &
& moment_integration(poly_order,tau,poly_multiplication(vpar_coef,x_vpar_coef)))
!
!print*, 'vpar:',vpar_coef, tau, moment_integration(poly_order,tau,vpar_coef)
!
optional_quantities%vpar_int = optional_quantities%vpar_int + dble(sign_rhs) * &
& moment_integration(poly_order,tau,poly_multiplication(dtdtau_coef,vpar_coef))
endif
!
!Optional computation of integral of v_par**2
if (boole_vpar2_int) then
! optional_quantities%vpar2_int = optional_quantities%vpar2_int + dt_dtau_const*&
! & moment_integration(poly_order,tau,poly_multiplication(vpar_coef,vpar_coef))
optional_quantities%vpar2_int = optional_quantities%vpar2_int + &
!
!First term
& hamiltonian_time(ind_tetr)%h1_in_curlA * &
& moment_integration(poly_order,tau,poly_multiplication(vpar_coef,vpar_coef)) + &
!
!Second term
& cm_over_e * hamiltonian_time(ind_tetr)%h1_in_curlh * &
& moment_integration(poly_order,tau,poly_multiplication(vpar_coef,poly_multiplication(vpar_coef,vpar_coef)))+ &
!
!Third term
& sum( hamiltonian_time(ind_tetr)%vec_mismatch_der * &
& moment_integration(poly_order,tau,poly_multiplication(poly_multiplication(vpar_coef,vpar_coef),x_coef)) ) + &
!
!Fourth term
& cm_over_e * sum(hamiltonian_time(ind_tetr)%vec_parcurr_der * &
& moment_integration(poly_order,tau,poly_multiplication(poly_multiplication(vpar_coef,vpar_coef),x_vpar_coef)))
!

!print*, 'vpar2', poly_multiplication(vpar_coef,vpar_coef), tau, &
!moment_integration(poly_order,tau,poly_multiplication(vpar_coef,vpar_coef))
!
optional_quantities%vpar2_int = optional_quantities%vpar2_int + dble(sign_rhs) * &
& moment_integration(poly_order,tau,poly_multiplication(dtdtau_coef,poly_multiplication(vpar_coef,vpar_coef)))
endif
!
deallocate(vpar_coef)
if(boole_time_hamiltonian.or.boole_vpar_int.or.boole_vpar2_int) deallocate(x_coef,x_vpar_coef)
deallocate(vpar_coef,x_coef,x_vpar_coef)
!
end subroutine calc_optional_quantities
!
Expand Down
2 changes: 2 additions & 0 deletions SRC/strong_electric_field_mod.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -99,6 +99,8 @@ subroutine get_electric_potential(R,phi,Z,electric_potential)
double precision :: d1,d2,d3,d4,d5,d6,d7,d8,d9,d10,d11,d12 !Dummy variables for field() to set psif
!
select case(electric_potential_option)
!
! create case (1) here where I can import data from Markus
!
!Get scalar electric potential as a rescale of the poloidal magnetic flux (potential contours align with flux surfaces)
case(2)
Expand Down
Loading