diff --git a/plasmaExpansion.f90 b/plasmaExpansion.f90 index 513d2d1..e15f053 100644 --- a/plasmaExpansion.f90 +++ b/plasmaExpansion.f90 @@ -354,13 +354,13 @@ program plasmaExpansion ! Set input parameters (remember these have to be in non-dimensional units) Temp0 = 60.0_dp * eV_to_K / Temp_ref - TempF = 5.0_dp * eV_to_K / Temp_ref + TempF = 60.0_dp * eV_to_K / Temp_ref n_ecr = 1.0e19_dp * cm3_to_m3 / n_ref c_s = sqrt(T_to_Z(Temp0) * gam * Temp0) u_bc0 = c_s!sqrt(Temp0) - u_bcF = sqrt(TempF) + u_bcF = u_bc0!sqrt(TempF) n_bc0 = n_ecr / T_to_Z(Temp0) - n_bcF = n_ecr*1.0e-1 / T_to_Z(Temp0) + n_bcF = n_bc0!n_ecr*1.0e-1 / T_to_Z(Temp0) ! Set domain boundaries (non-dimensional units) r0 = 200.0e-6_dp / L_ref @@ -395,7 +395,7 @@ program plasmaExpansion t0 = 0.0_dp tf = 1.0e-6_dp / t_ref ! tf = 1.0e1_dp * (rf - r0) / c_s - dt = 1.0e-2_dp*dr/c_s + dt = 1.0e-2_dp*dr/(10.0_dp*c_s) nt = nint((tf - t0) / dt) dt = (tf - t0) / float(nt) @@ -430,7 +430,7 @@ program plasmaExpansion u_i = 0.0_dp E_i = 0.0_dp T_i = 0.0_dp - n_e = 0.0_dp + n_e = 1.0_dp Zave = 0.0_dp phi = 0.0_dp phi_old = 0.0_dp @@ -445,14 +445,14 @@ program plasmaExpansion diag_high = 0.0_dp b = 0.0_dp diag = -2.0_dp / dr**2 - diag_low = 1.0_dp / dr**2 * r(1:nr-1) / r(2:nr) - diag_high = 1.0_dp / dr**2 * r(2:nr) / r(1:nr-1) - ! diag(1) = 1.0_dp ! Dirichlet - ! diag_high(1) = 0.0_dp ! Dirichlet - diag_high(1) = 2.0_dp / dr**2 ! Neumann - diag(nr) = 1.0_dp ! Dirichlet - diag_low(nr-1) = 0.0_dp ! Dirichlet - ! diag_low(nr-1) = 2.0_dp / dr**2 ! Neumann + diag_low = 1.0_dp / dr**2! * r(1:nr-1) / r(2:nr) + diag_high = 1.0_dp / dr**2! * r(2:nr) / r(1:nr-1) + diag(1) = 1.0_dp ! Dirichlet + diag_high(1) = 0.0_dp ! Dirichlet + ! diag_high(1) = 2.0_dp / dr**2 ! Neumann + ! diag(nr) = 1.0_dp ! Dirichlet + ! diag_low(nr-1) = 0.0_dp ! Dirichlet + diag_low(nr-1) = 2.0_dp / dr**2 ! Neumann allocate(A(1:nr,1:nr)) A = 0.0_dp @@ -470,10 +470,9 @@ program plasmaExpansion Res = 0.0_dp ! Set boundary values - ! phi0 = 0.0_dp / phi_ref ! Dirichlet - ! phi(1) = phi0 ! Dirichlet - phi0 = phi(1) ! Neumann - phiF = 0.0_dp ! Dirichlet + phi0 = 0.0_dp / phi_ref ! Dirichlet + ! phi0 = phi(1) ! Neumann + ! phiF = 0.0_dp ! Dirichlet allocate(f0(j0:nv)) f0 = 0.0_dp @@ -520,13 +519,13 @@ program plasmaExpansion 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)) + f_i(i,1:j0-1) = f_i_old(i,1:j0-1) - v(1:j0-1)*dt/dr*(f_i_old(i+1,1:j0-1) - & + 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)) + f_i(i,j0:nv) = f_i_old(i, j0:nv) - v( j0:nv)*dt/dr*(f_i_old(i , j0:nv) - & + f_i_old(i-1, j0:nv)) end if n_i(i) = sum(f_i(i,:)*dv) @@ -546,29 +545,26 @@ program plasmaExpansion !$omp end parallel do ! Solve Poission (maximum number of iterations, break if convergence is reached before) - do k = 1, 500 + do k = 1, 1000 ! Store previous value phi_old = phi - ! Calculate distribution of electrons - n_e = Zave(1) * n_i(1) * exp((phi_old - phi0) / T_i(1)) - ! Diagonal matrix for Newton integration scheme diag = -2.0_dp / dr**2 - n_e - diag_low = 1.0_dp / dr**2 * r(1:nr-1) / r(2:nr) - diag_high = 1.0_dp / dr**2 * r(2:nr) / r(1:nr-1) - ! diag(1) = 1.0_dp ! Dirichlet - ! diag_high(1) = 0.0_dp ! Dirichlet - diag_high(1) = 2.0_dp / dr**2 - n_e(1) ! Neumann - diag(nr) = 1.0_dp ! Dirichlet - diag_low(nr-1) = 0.0_dp ! Dirichlet - ! diag_low(nr-1) = 2.0_dp / dr**2 - n_e(nr) ! Neumann + diag_low = 1.0_dp / dr**2! * r(1:nr-1) / r(2:nr) + diag_high = 1.0_dp / dr**2! * r(2:nr) / r(1:nr-1) + diag(1) = 1.0_dp ! Dirichlet + diag_high(1) = 0.0_dp ! Dirichlet + ! diag_high(1) = 2.0_dp / dr**2 - n_e(1) ! Neumann + ! diag(nr) = 1.0_dp ! Dirichlet + ! diag_low(nr-1) = 0.0_dp ! Dirichlet + diag_low(nr-1) = 2.0_dp / dr**2 - n_e(nr) ! Neumann ! Calculate charge density b = -(Zave*n_i - n_e) ! Apply boundary conditions - ! b(1) = phi0 ! Dirichlet - b(nr) = phiF ! Dirichlet + b(1) = phi0 ! Dirichlet + ! b(nr) = phiF ! Dirichlet ! Calculate residual !$omp parallel workshare @@ -578,39 +574,39 @@ program plasmaExpansion ! Iterate system call dgtsv(nr, 1, diag_low, diag, diag_high, Res, nr, info) phi = phi_old + Res - phi(1) = phi(2) ! Ensures smooth transition in Neumann boundary - phi0 = phi(1) ! Neumann + + ! Calculate distribution of electrons + n_e = Zave(1) * n_i(1) * exp((phi - phi0) / T_i(1)) ! Check if the solution has converged phiConv = maxval(abs(Res),1) - if (phiConv < 1.0e-1_dp) then + if (phiConv < 1.0e-5_dp) then exit end if - ! Calculate new potential to ensure 0 current at the edge - if (n_i(nr) > 1.0e-10_dp) then - phiF = phi0 + T_i(1) * log((2.0_dp*sqrt(pi)*Zave(nr)*n_i(nr)*u_i(nr)) / (Zave(1)*n_i(1)*sqrt(m_i*T_i(1)/m_e))) - - else - phiF = phi(nr-5) - - end if + ! ! Calculate new potential to ensure 0 current at the edge + ! if (n_i(nr) > 1.0e-10_dp) then + ! phiF = phi0 + T_i(1) * log((2.0_dp*sqrt(pi)*Zave(nr)*n_i(nr)*u_i(nr)) / (Zave(1)*n_i(1)*sqrt(m_i*T_i(1)/m_e))) + ! + ! else + ! phiF = phi(nr-5) + ! + ! end if end do ! Calculate electric field - ! E(1) = - (phi(2) - phi(1)) / dr ! Dirichlet - E(1) = 0.0_dp ! Neumann + E(1) = - (phi(2) - phi(1)) / dr ! Dirichlet + ! E(1) = 0.0_dp ! Neumann !$omp parallel do do i = 2, nr-1 E(i) = - 0.5_dp*(phi(i+1) - phi(i-1)) / dr end do !$omp end parallel do - E(nr) = - (phi(nr) - phi(nr-1)) / dr ! Dirichlet - ! E(nr) = 0.0_dp ! Neumann - ! Trick to avoid problems at the sheath + ! E(nr) = - (phi(nr) - phi(nr-1)) / dr ! Dirichlet + E(nr) = 0.0_dp ! Neumann ! Update intermediate f f_i_old = f_i