Start to implement options for collision operator and COSOlver loader

src/array_mod.F90

@@ -8,6 +8,10 @@ MODULE array
   COMPLEX(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: moments_rhs_e ! (ip,ij,ikr,ikz,updatetlevel)
   COMPLEX(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: moments_rhs_i ! (ip,ij,ikr,ikz,updatetlevel)
 
+  ! To load collision matrix
+  REAL(dp), DIMENSION(:,:), ALLOCATABLE :: iCe ! (ip,ij)
+  REAL(dp), DIMENSION(:,:), ALLOCATABLE :: iCi ! (ip,ij)
+
   ! Intermediate steps in rhs of equations
   !COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE:: moments_Apl, moments_Bpl, moments_Cpl, moments_Dpl,&
   !                                        moments_Epl, moments_Fpl, moments_Gpl, moments_Hpl  

src/inital.F90

 SUBROUTINE inital
   
-  use basic
-
-  use prec_const
+  USE basic
+  USE model, ONLY : CO
+  USE prec_const
   implicit none
 
   WRITE(*,'(a/)') '=== Set initial conditions ==='
 
   CALL init_profiles
 
+  IF (CO .EQ. -1) THEN
+    WRITE(*,'(a/)') '=== Load Full Coulomb matrix ==='
+
+    CALL load_FC_mat
+  ENDIF
+  !
+
 END SUBROUTINE inital
 
 
@@ -18,10 +25,10 @@ SUBROUTINE init_profiles
   USE basic
   USE fourier_grid
   USE fields
-  use initial_par
+  USE initial_par
   USE time_integration
 
-  use prec_const
+  USE prec_const
   IMPLICIT NONE
 
   INTEGER :: ip,ij, ikr,ikz
@@ -62,3 +69,9 @@ SUBROUTINE init_profiles
 
 END SUBROUTINE init_profiles
 
+SUBROUTINE load_FC_mat
+  USE basic
+  USE array
+  IMPLICIT NONE
+
+END SUBROUTINE load_FC_mat
\ No newline at end of file

src/memory.F90

@@ -6,8 +6,9 @@ SUBROUTINE memory
   USE fields
   USE fourier_grid
   USE time_integration  
+  USE model, ONLY: CO
 
-  use prec_const
+  USE prec_const
   IMPLICIT NONE
 
   ! Moments and moments rhs
@@ -19,4 +20,9 @@ SUBROUTINE memory
   ! Electrostatic potential
   CALL allocate_array(phi, ikrs,ikre, ikzs,ikze)
 
+  ! Collision matrix
+  IF (CO .EQ. -1) THEN
+    CALL allocate_array( iCe, ips_i,ipe_i, ijs_i,ije_i )
+    CALL allocate_array( iCi, ips_i,ipe_i, ijs_i,ije_i )
+  ENDIF
 END SUBROUTINE memory

src/model_mod.F90

@@ -5,17 +5,18 @@ MODULE model
   IMPLICIT NONE
   PRIVATE
 
-  REAL(dp), PUBLIC, PROTECTED ::      nu = 1._dp     ! Collision frequency
-  REAL(dp), PUBLIC, PROTECTED ::   tau_e = 1._dp     ! Collision frequency
-  REAL(dp), PUBLIC, PROTECTED ::   tau_i = 1._dp     ! Collision frequency
-  REAL(dp), PUBLIC, PROTECTED :: sigma_e = 1._dp     ! Collision frequency
-  REAL(dp), PUBLIC, PROTECTED :: sigma_i = 1._dp     ! Collision frequency
-  REAL(dp), PUBLIC, PROTECTED ::     q_e = 1._dp     ! Collision frequency
-  REAL(dp), PUBLIC, PROTECTED ::     q_i = 1._dp     ! Collision frequency
-  REAL(dp), PUBLIC, PROTECTED ::   eta_n = 1._dp     ! Collision frequency
-  REAL(dp), PUBLIC, PROTECTED ::   eta_T = 1._dp     ! Collision frequency
-  REAL(dp), PUBLIC, PROTECTED ::   eta_B = 1._dp     ! Collision frequency
-  REAL(dp), PUBLIC, PROTECTED :: lambdaD = 1._dp     ! Collision frequency
+  INTEGER(dp), PUBLIC, PROTECTED ::   CO =  0         ! Collision Operator
+  REAL(dp), PUBLIC, PROTECTED ::      nu =  1._dp     ! Collision frequency
+  REAL(dp), PUBLIC, PROTECTED ::   tau_e =  1._dp     ! Temperature
+  REAL(dp), PUBLIC, PROTECTED ::   tau_i =  1._dp     !
+  REAL(dp), PUBLIC, PROTECTED :: sigma_e =  1._dp     ! Mass
+  REAL(dp), PUBLIC, PROTECTED :: sigma_i =  1._dp     !
+  REAL(dp), PUBLIC, PROTECTED ::     q_e = -1._dp     ! Charge
+  REAL(dp), PUBLIC, PROTECTED ::     q_i =  1._dp     !
+  REAL(dp), PUBLIC, PROTECTED ::   eta_n =  1._dp     ! Density gradient
+  REAL(dp), PUBLIC, PROTECTED ::   eta_T =  1._dp     ! Temperature gradient
+  REAL(dp), PUBLIC, PROTECTED ::   eta_B =  1._dp     ! Magnetic gradient
+  REAL(dp), PUBLIC, PROTECTED :: lambdaD =  1._dp     ! Debye length
 
   PUBLIC :: model_readinputs, model_outputinputs
 
@@ -28,7 +29,7 @@ CONTAINS
     USE prec_const
     IMPLICIT NONE
 
-    NAMELIST /MODEL_PAR/ nu, tau_e, tau_i, sigma_e, sigma_i, &
+    NAMELIST /MODEL_PAR/ CO, nu, tau_e, tau_i, sigma_e, sigma_i, &
                          q_e, q_i, eta_n, eta_T, eta_B, lambdaD
 
     READ(lu_in,model_par)
@@ -49,6 +50,7 @@ CONTAINS
 
     INTEGER, INTENT(in) :: fidres
     CHARACTER(len=256), INTENT(in) :: str
+    CALL attach(fidres, TRIM(str),      "CO",      CO)
     CALL attach(fidres, TRIM(str),      "nu",      nu)
     CALL attach(fidres, TRIM(str),   "tau_e",   tau_e)
     CALL attach(fidres, TRIM(str),   "tau_i",   tau_i)

wk/kz_linear.m

@@ -13,18 +13,19 @@ OUTPUTS.write_phi     = '.true.';
 OUTPUTS.write_doubleprecision = '.true.';
 OUTPUTS.resfile0      = '''results''';
 %% Grid parameters
-GRID.pmaxe = 81;
-GRID.jmaxe = 20;
-GRID.pmaxi = 81;
-GRID.jmaxi = 20;
+GRID.pmaxe = 25;
+GRID.jmaxe = 9;
+GRID.pmaxi = 25;
+GRID.jmaxi = 9;
 GRID.nkr   = 1;
 GRID.krmin = 0.;
 GRID.krmax = 0.;
-GRID.nkz   = 10;
+GRID.nkz   = 20;
 GRID.kzmin = 0.1/sqrt(2.0);
-GRID.kzmax = 2.8/sqrt(2.0);
+GRID.kzmax = 5.0/sqrt(2.0);
 %% Model parameters
-MODEL.nu      = 0.001;
+MODEL.CO      = -1;  % Collision operator (-1 = Full Coulomb, 0 = Dougherty)
+MODEL.nu      = 0.001; % collisionality nu*L_perp/Cs0
 MODEL.tau_e   = 1.0;
 MODEL.tau_i   = 1.0;
 MODEL.sigma_e = 0.0233380;