From 3754c0b910fde7c76a297be71831d94c1fb65609 Mon Sep 17 00:00:00 2001 From: JHendrikx Date: Mon, 3 Feb 2025 20:56:04 +0100 Subject: [PATCH] Add additional nZ rank to arrays --- vlaplex.f90 | 178 +++++++++++++++++++++++++++++----------------------- 1 file changed, 98 insertions(+), 80 deletions(-) diff --git a/vlaplex.f90 b/vlaplex.f90 index 6053586..d82a49d 100644 --- a/vlaplex.f90 +++ b/vlaplex.f90 @@ -46,8 +46,8 @@ program VlaPlEx real(dp):: t0, tf real(dp):: time real(dp):: dr, dv, dt - integer:: nr, nv, nt - integer:: i, j, t + integer:: nr, nv, nt, nz + integer:: i, iz, j, t integer:: j0 ! First integer of positive velocity real(dp):: Temp_bc ! Temperature @@ -58,12 +58,12 @@ program VlaPlEx type(tableBC):: boundaryConditions character(:), allocatable:: bc_file - real(dp), allocatable, dimension(:,:):: f_i, f_i_old - real(dp), allocatable, dimension(:):: f0 ! Boundary at r = x_0 - real(dp), allocatable, dimension(:):: n_i - real(dp), allocatable, dimension(:):: u_i + real(dp), allocatable, dimension(:,:,:):: f_i, f_i_old + real(dp), allocatable, dimension(:,:):: f0 ! Boundary at r = x_0 + real(dp), allocatable, dimension(:,:):: n_i + real(dp), allocatable, dimension(:,:):: u_i real(dp), allocatable, dimension(:):: E_i - real(dp), allocatable, dimension(:):: T_i + real(dp), allocatable, dimension(:,:):: T_i real(dp), allocatable, dimension(:):: n_e real(dp), allocatable, dimension(:):: Zave real(dp), allocatable, dimension(:):: diag, diag_low, diag_high @@ -79,7 +79,7 @@ program VlaPlEx ! real(dp):: phiF integer:: k - real(dp), allocatable, dimension(:):: fCum_i + real(dp), allocatable, dimension(:,:):: fCum_i real(dp):: rCum integer:: rCum_index @@ -96,7 +96,7 @@ program VlaPlEx ! Set input parameters (remember these have to be in non-dimensional units) c_s = sqrt(11.0_dp * gamma_i * 1.0_dp) - bc_file = 'bc_80ns_T60Z16.csv' + bc_file = 'bc.csv' call boundaryConditions%init(bc_file) ! Set domain boundaries (non-dimensional units) @@ -156,14 +156,15 @@ program VlaPlEx write(*, '(A,ES0.4e3)') 'CFL: ', dt*vf/dr + nz = 2 ! Allocate vectors - allocate(f_i(1:nr,1:nv), f_i_old(1:nr,1:nv)) - allocate(n_i(1:nr)) - allocate(u_i(1:nr), E_i(1:nr), T_i(1:nr)) + allocate(f_i(1:nr,1:nv,1:nz), f_i_old(1:nr,1:nv,1:nz)) + allocate(n_i(1:nr,1:nz)) + allocate(u_i(1:nr,1:nz), E_i(1:nr), T_i(1:nr,1:nz)) allocate(Zave(1:nr)) allocate(n_e(1:nr)) allocate(phi(1:nr), phi_old(1:nr), E(1:nr)) - allocate(fCum_i(1:nv)) + allocate(fCum_i(1:nv,1:nz)) f_i = 0.0_dp f_i_old = 0.0_dp n_i = 0.0_dp @@ -216,7 +217,7 @@ program VlaPlEx phi0 = 1.0e2_dp / phi_ref ! Dirichlet phi(1) = phi0 ! Dirichlet ! phi0 = phi(1) ! Neumann - allocate(f0(j0:nv)) + allocate(f0(j0:nv,1:nz)) f0 = 0.0_dp ! Output initial values @@ -227,7 +228,7 @@ program VlaPlEx ! call writeOutputF(t, dt, nr, r, nv, v, f_i_old) call writeOutputFCum(t, dt, r(rCum_index), nv, v, fCum_i) call writeOutputPhi(t, dt, nr, r, phi, E, n_e) - call writeOutputMom(t, dt, nr, r, n_i, u_i, T_i, Zave) + call writeOutputMom(t, dt, nr, r, n_i(:,1), u_i(:,1), T_i(:,1), Zave) ! Main loop do t = 1, nt @@ -235,9 +236,17 @@ program VlaPlEx call boundaryConditions%get(time, n_bc, u_bc, Temp_bc, Zave_bc) call writeOutputBoundary(t, dt, n_bc, u_bc, Temp_bc, Zave_bc) u_bc = sqrt(Zave_bc * Temp_bc) - ! f0(j0:nv) = v(j0:nv)**2 / sqrt(PI*Temp_bc**3) * exp(-(v(j0:nv) - u_bc)**2 / Temp_bc) - f0(j0:nv) = 1.0_dp / sqrt(PI*Temp_bc) * exp(-(v(j0:nv) - u_bc)**2 / Temp_bc) - f0 = f0 * n_bc / (sum(f0)*dv) + do iz = 1, nz + ! f0(j0:nv) = v(j0:nv)**2 / sqrt(PI*Temp_bc**3) * exp(-(v(j0:nv) - u_bc)**2 / Temp_bc) + f0(j0:nv,iz) = 1.0_dp / sqrt(PI*Temp_bc) * exp(-(v(j0:nv) - u_bc)**2 / Temp_bc) + f0(:,iz) = f0(:,iz) * n_bc / (sum(f0(:,iz))*dv) + + ! Boundary conditions + ! r = r0, v>0 + f_i_old(1,j0:nv,iz) = f0(:,iz) + f_i(1,j0:nv,iz) = f_i_old(1,j0:nv,iz) + T_i(1,iz) = Temp_bc + end do T_e = Temp_bc print *, 'Time: ', time * t_ref print *, 'Temp_bc: ', Temp_bc @@ -246,50 +255,51 @@ program VlaPlEx print *, '-------------------------' - ! Boundary conditions - ! r = r0, v>0 - f_i_old(1,j0:nv) = f0 - f_i(1,j0:nv) = f_i_old(1,j0:nv) - T_i(1) = Temp_bc + Zave(1) = Zave_bc ! r = rf, v<0 - f_i_old(nr,1:j0-1) = 0.0_dp - f_i(nr,1:j0-1) = f_i_old(nr,1:j0-1) + f_i_old(nr,1:j0-1,:) = 0.0_dp + f_i(nr,1:j0-1,:) = f_i_old(nr,1:j0-1,:) ! set edge velocities to 0 - f_i_old(:,1) = 0.0_dp - f_i_old(:,nv) = 0.0_dp + f_i_old(:,1,:) = 0.0_dp + f_i_old(:,nv,:) = 0.0_dp ! Advect in the r direction !$omp parallel do - do i = 1, nr - ! Advect negative velocity - if (i < nr) then - f_i(i,1:j0-1) = f_i_old(i,1:j0-1) - v(1:j0-1)*dt/dr/r(i)**2*(r(i+1)**2*f_i_old(i+1,1:j0-1) - & - r(i )**2*f_i_old(i ,1:j0-1)) - end if - ! Advect positive velocity - if (i > 1) then - f_i(i,j0:nv) = f_i_old(i, j0:nv) - v( j0:nv)*dt/dr/r(i)**2*(r(i )**2*f_i_old(i , j0:nv) - & - r(i-1)**2*f_i_old(i-1, j0:nv)) - end if + do iz = 1, nz + do i = 1, nr + ! Advect negative velocity + if (i < nr) then + f_i(i,1:j0-1,iz) = f_i_old(i,1:j0-1,iz) - v(1:j0-1)*dt/dr/r(i)**2*(r(i+1)**2*f_i_old(i+1,1:j0-1,iz) - & + r(i )**2*f_i_old(i ,1:j0-1,iz)) + end if + ! Advect positive velocity + if (i > 1) then + f_i(i,j0:nv,iz) = f_i_old(i, j0:nv,iz) - v( j0:nv)*dt/dr/r(i)**2*(r(i )**2*f_i_old(i , j0:nv,iz) - & + r(i-1)**2*f_i_old(i-1, j0:nv,iz)) + end if - n_i(i) = sum(f_i(i,:))*dv - if (n_i(i) > 1.0e-10_dp) then - u_i(i) = sum(v(:) *f_i(i,:))*dv / n_i(i) - E_i(i) = sum(v(:)**2*f_i(i,:))*dv / n_i(i) - T_i(i) = 2.0_dp*E_i(i) - 2.0_dp*u_i(i)**2 - Zave(i) = Zave_bc + n_i(i,iz) = sum(f_i(i,:,iz))*dv + if (n_i(i,1) > 1.0e-10_dp) then + u_i(i,iz) = sum(v(:) *f_i(i,:,iz))*dv / n_i(i,iz) + E_i(i) = sum(v(:)**2*f_i(i,:,iz))*dv / n_i(i,iz) + T_i(i,iz) = 2.0_dp*E_i(i) - 2.0_dp*u_i(i,iz)**2 + Zave(i) = Zave_bc - else - u_i(i) = 0.0_dp - T_i(i) = 0.0_dp - Zave(i) = 0.0_dp - end if + else + u_i(i,iz) = 0.0_dp + T_i(i,iz) = 0.0_dp + Zave(i) = 0.0_dp + end if + end do end do !$omp end parallel do ! Assume quasi-neutrality to start iterating - n_e = Zave * n_i + n_e = Zave * n_i(:,1) + !do iz = 1, nz + ! n_e = n_e + Zave * n_i(:,iz) + !end do db_dphi = 0.0_dp ! Solve Poission (maximum number of iterations, break if convergence is reached before) @@ -307,7 +317,12 @@ program VlaPlEx diag_low(nr-1) = 2.0_dp / dr**2 - db_dphi(nr) ! Neumann ! Calculate charge density - b = -(Zave*n_i - n_e) + !b = n_e + !b = b - (Zave * n_i(:,1)) + !do iz = 1, nz + ! b = b - (Zave * n_i(:,iz)) + !end do + b = -(Zave * n_i(:,1) - n_e) ! Apply boundary conditions b(1) = phi0 ! Dirichlet ! b(nr) = 0.0_dp ! Dirichlet @@ -324,10 +339,10 @@ program VlaPlEx ! Calculate distribution of electrons ! n_e = Zave(1) * n_i(1) * exp((phi- phi0) / T_e) ! Isothermal (Boltzmann) - n_e = Zave(1) * n_i(1) * (1.0_dp + (gamma_e - 1.0_dp)/gamma_e*(phi-phi0)/T_e)**gamma_e_exp !Polytropic + n_e = Zave(1) * n_i(1,1) * (1.0_dp + (gamma_e - 1.0_dp)/gamma_e*(phi-phi0)/T_e)**gamma_e_exp !Polytropic ! Diagonal matrix for Newton integration scheme ! db_dphi = n_e / T_e ! Isotropic - db_dphi = Zave(1) * n_i(1) / (gamma_e * T_e) * & + db_dphi = Zave(1) * n_i(1,1) / (gamma_e * T_e) * & (1.0_dp + (gamma_e - 1.0_dp)/gamma_e*(phi-phi0)/T_e)**gamma_e_dexp !Polytropic ! Check if the solution has converged @@ -363,46 +378,49 @@ program VlaPlEx ! Update intermediate f f_i_old = f_i - ! Advect in the v direction - ! i = 1, v<0 - i = 1 - if (E(i) >= 0.0_dp) then - f_i(i,2:j0-2) = f_i_old(i,2:j0-2) - Zave(i)*E(i)*dt/dv*(f_i_old(i,2:j0-2) - f_i_old(i,1:j0-3)) - else - f_i(i,2:j0-2) = f_i_old(i,2:j0-2) - Zave(i)*E(i)*dt/dv*(f_i_old(i,3:j0-1) - f_i_old(i,2:j0-2)) - - end if - ! i = 2, nr-1; all v - !$omp parallel do - do i = 2, nr-1 + do iz = 1, nz + ! Advect in the v direction + ! i = 1, v<0 + i = 1 if (E(i) >= 0.0_dp) then - f_i(i,2:nv-1) = f_i_old(i,2:nv-1) - Zave(i)*E(i)*dt/dv*(f_i_old(i,2:nv-1) - f_i_old(i,1:nv-2)) + f_i(i,2:j0-2,iz) = f_i_old(i,2:j0-2,iz) - Zave(i)*E(i)*dt/dv*(f_i_old(i,2:j0-2,iz) - f_i_old(i,1:j0-3,iz)) else - f_i(i,2:nv-1) = f_i_old(i,2:nv-1) - Zave(i)*E(i)*dt/dv*(f_i_old(i,3:nv) - f_i_old(i,2:nv-1)) + f_i(i,2:j0-2,iz) = f_i_old(i,2:j0-2,iz) - Zave(i)*E(i)*dt/dv*(f_i_old(i,3:j0-1,iz) - f_i_old(i,2:j0-2,iz)) end if + ! i = 2, nr-1; all v + !$omp parallel do + do i = 2, nr-1 + if (E(i) >= 0.0_dp) then + f_i(i,2:nv-1,iz) = f_i_old(i,2:nv-1,iz) - Zave(i)*E(i)*dt/dv*(f_i_old(i,2:nv-1,iz) - f_i_old(i,1:nv-2,iz)) + else + f_i(i,2:nv-1,iz) = f_i_old(i,2:nv-1,iz) - Zave(i)*E(i)*dt/dv*(f_i_old(i,3:nv,iz) - f_i_old(i,2:nv-1,iz)) + end if + + end do + !$omp end parallel do + ! i = nr, v>=0 + i = nr + if (E(i) >= 0.0_dp) then + f_i(i,j0+1:nv-1,iz) = f_i_old(i,j0+1:nv-1,iz) - Zave(i)*E(i)*dt/dv*(f_i_old(i,j0+1:nv-1,iz) - f_i_old(i,j0:nv-2,iz)) + else + f_i(i,j0+1:nv-1,iz) = f_i_old(i,j0+1:nv-1,iz) - Zave(i)*E(i)*dt/dv*(f_i_old(i,j0+2:nv,iz) - f_i_old(i,j0+1:nv-1,iz)) + + end if end do - !$omp end parallel do - ! i = nr, v>=0 - i = nr - if (E(i) >= 0.0_dp) then - f_i(i,j0+1:nv-1) = f_i_old(i,j0+1:nv-1) - Zave(i)*E(i)*dt/dv*(f_i_old(i,j0+1:nv-1) - f_i_old(i,j0:nv-2)) - else - f_i(i,j0+1:nv-1) = f_i_old(i,j0+1:nv-1) - Zave(i)*E(i)*dt/dv*(f_i_old(i,j0+2:nv) - f_i_old(i,j0+1:nv-1)) - - end if ! Reset values for next iteration f_i_old = f_i - fCum_i = fCum_i + f_i_old(rCum_index,:) - + do iz = 1, nz + fCum_i(:,iz) = fCum_i(:,iz) + f_i_old(rCum_index,:,iz) + end do ! Write output if (mod(t,everyOutput) == 0 .or. t == nt) then ! call writeOutputF(t, dt, nr, r, nv, v, f_i_old) call writeOutputPhi(t, dt, nr, r, phi, E, n_e) - call writeOutputMom(t, dt, nr, r, n_i, u_i, T_i, Zave) - call writeOutputFCum(t, dt, r(rCum_index), nv, v, fCum_i) + call writeOutputMom(t, dt, nr, r, n_i(:,1), u_i(:,1), T_i(:,1), Zave) + call writeOutputFCum(t, dt, r(rCum_index), nv, v, fCum_i(:,1)) end if