Implemented the 2nd interpolation for species and field

.gitignore

@@ -10,3 +10,4 @@
 .depend
 bin
 source/build
+docs/build

config/make.gnu_openmpi

@@ -1,6 +1,6 @@
 # Fortran compiler configuration
 FC = mpif90
-FC_OPTS = -g -O3 -fdefault-real-8 -fdefault-double-8 -fopenmp -ffree-form 
+FC_OPTS = -Wa,-q -cpp -g -O3 -fdefault-real-8 -fdefault-double-8 -fopenmp -ffree-form 
 # -fbacktrace -fbounds-check
 # -fbounds-check
 # FC_OPTS = -g -Wall -Wextra -Wconversion -fimplicit-none -fbacktrace \
@@ -12,7 +12,7 @@ CC = gcc
 CC_OPTS = -O3 -std=c99
 
 # Preprocessor configuration
-FPP = gcc -C -P -E -cpp
+FPP = gcc -C -P -E -x assembler-with-cpp
 
 # Python 2.x configuration. Used for automatically generate compile rules
 PYTHON = python 

config/make.template.gnu

@@ -10,7 +10,7 @@
 FC = mpif90
 
 # Configure the compiler options
-FC_OPTS = -g -O3 -fdefault-real-8 -fdefault-double-8 -fopenmp -ffree-form
+FC_OPTS = -Wa,-q -cpp -g -O3 -fdefault-real-8 -fdefault-double-8 -fopenmp -ffree-form
 
 # Enable the generation of run-time checks. The argument shall be a comma-delimited 
 # list of the following keywords:
@@ -45,7 +45,7 @@ CC = gcc
 CC_OPTS = -O3 -std=c99
 
 # Preprocessor configuration
-FPP = gcc -C -P -E -cpp
+FPP = gcc -C -P -E -x assembler-with-cpp 
 
 # ------------------------------------------------------------------------------
 # Configure the modules that will NOT be installed

source/Makefile

@@ -90,11 +90,13 @@ $(main_obj) : $(configfile) | $(builddir)
 $(builddir) :
 	@echo "Creating $@ directory"
 	@mkdir -p $@
+	@mkdir -p $@/preprocessed
 
 $(builddir)/%.o : %.f03
 	@echo "[F03] $(<F)"
 	@$(FPP) $(FPPF) -I. -I./$(builddir) $< > $(builddir)/$(<F)
-	@cd $(builddir) && $(FC) $(INCF) -c $(FC_OPTS) $(<F) -o $(@F) && rm $(<F)
+	@cd $(builddir) && $(FC) $(INCF) -c $(FC_OPTS) $(<F) 
+	@cd $(builddir) && mv $(<F) ./preprocessed
 
 $(builddir)/%.o : %.c
 	@echo "[CC] $(<F)"

source/fields/field_b_class.f03

@@ -112,7 +112,9 @@ subroutine init_field_b( this, opts, max_mode, part_shape, boundary, relax_fac,
 
   case ( p_ps_quadratic )
 
-    call write_err( "Quadratic particle shape not implemented." )
+    gc_num(:,1) = (/1, 2/)
+    gc_num(:,2) = (/0, 1/)
+    ! call write_err( "Quadratic particle shape not implemented." )
 
   case default
 

source/fields/field_e_class.f03

@@ -88,7 +88,9 @@ subroutine init_field_e( this, opts, max_mode, part_shape, boundary, entity )
 
   case ( p_ps_quadratic )
 
-    call write_err( "Quadratic particle shape not implemented." )
+    gc_num(:,1) = (/1, 2/)
+    gc_num(:,2) = (/0, 1/)
+    ! call write_err( "Quadratic particle shape not implemented." )
 
   case default
 

source/fields/field_psi_class.f03

@@ -79,7 +79,10 @@ subroutine init_field_psi( this, opts, max_mode, part_shape, boundary )
 
   case ( p_ps_quadratic )
 
-    call write_err( "Quadratic particle shape not implemented." )
+
+    gc_num(:,1) = (/1, 2/)
+    gc_num(:,2) = (/0, 1/)
+    ! call write_err( "Quadratic particle shape not implemented." )
 
   case default
 

source/fields/field_src_class.f03

@@ -81,8 +81,11 @@ subroutine init_field_rho(this, opts, max_mode, part_shape, smooth_order, has_2d
       gc_num(:,2) = (/0, 1/)
 
     case ( p_ps_quadratic )
-      call write_err( "Quadratic particle shape not implemented." )
-
+
+    !   call write_err( "Quadratic particle shape not implemented." )
+        gc_num(:,1) = (/1, 2/)
+        gc_num(:,2) = (/0, 2/)
+
     case default
       call write_err( "Invalid particle shape." )
 
@@ -149,8 +152,11 @@ subroutine init_field_jay(this, opts, max_mode, part_shape, smooth_order)
       gc_num(:,2) = (/0, 1/)
 
     case ( p_ps_quadratic )
-
-      call write_err( "Quadratic particle shape not implemented." )
+
+      gc_num(:,1) = (/1, 2/)
+      gc_num(:,2) = (/0, 2/)
+
+    !   call write_err( "Quadratic particle shape not implemented." )
 
     case default
 
@@ -223,8 +229,11 @@ subroutine init_field_djdxi( this, opts, max_mode, part_shape, smooth_order )
     gc_num(:,2) = (/0, 1/)
 
   case ( p_ps_quadratic )
-
-    call write_err( "Quadratic particle shape not implemented." )
+
+    gc_num(:,1) = (/1, 2/)
+    gc_num(:,2) = (/0, 1/)
+
+    ! call write_err( "Quadratic particle shape not implemented." )
 
   case default
 

source/fields/field_vpot_class.f03

@@ -97,7 +97,9 @@ subroutine init_field_vpot( this, opts, max_mode, part_shape, boundary )
 
   case ( p_ps_quadratic )
 
-    call write_err( "Quadratic particle shape not implemented." )
+    gc_num(:,1) = (/1, 2/)
+    gc_num(:,2) = (/0, 1/)
+    ! call write_err( "Quadratic particle shape not implemented." )
 
   case default
 

source/interpolation.f03

@@ -29,14 +29,20 @@ subroutine spline_quadratic( x, s )
   real, intent(in) :: x
   real, dimension(-1:1), intent(inout) :: s
 
-  real :: t0, t1
+! Previous work, unfinished
+!   real :: t0, t1
 
-  t0 = 1.0 - x
-  t1 = x
+!   t0 = 1.0 - x
+!   t1 = x
 
-  s(-1) = 0.5 * t0**2
-  s( 0) = 0.5 + t0*t1
-  s( 1) = 0.5 * t1**2
+!   s(-1) = 0.5 * t0**2
+!   s( 0) = 0.5 + t0*t1
+!   s( 1) = 0.5 * t1**2
+
+  ! the x is in [-0.5, 0.5)
+  s(-1) = 0.5*(0.5-x)**2
+  s( 0) = 3.0/4.0 - x**2
+  s( 1) = 0.5*(0.5+x)**2
 
 end subroutine spline_quadratic
 

source/sim_beams_class.f03

@@ -95,8 +95,11 @@ subroutine init_sim_beams( this, input, opts )
   select case ( trim(str) )
   case ( 'linear' )
     ps = p_ps_linear
+  case ( 'quadratic' )
+    ps = p_ps_quadratic
+    call write_stdout('Attention! Quadratic interpolation for beams not implemented!')
   case default
-    call write_err( 'Invalid interpolation type! Only "linear" are supported &
+    call write_err( 'Invalid interpolation type! Only "linear" and "quadratic" are supported &
       &currently.' )
   end select
 

source/sim_fields_class.f03

@@ -128,8 +128,11 @@ subroutine init_sim_fields( this, input, opts )
   select case ( trim(str) )
   case ( 'linear' )
     ps = p_ps_linear
+  case ( 'quadratic' )
+    ps = p_ps_quadratic
+    ! call write_stdout('Attention! Quadratic interpolation for fields not implemented!')
   case default
-    call write_err( 'Invalid interpolation type! Only "linear" are supported currently.' )
+    call write_err( 'Invalid interpolation type! Only "linear" and "quadratic" are supported currently.' )
   end select
 
   relax_fac = 1.0d-3

source/sim_lasers_class.f03

@@ -91,6 +91,7 @@ subroutine init_sim_lasers( this, input, opts )
   ! local data
   character(len=32), save :: sname = 'init_sim_lasers'
   character(len=:), allocatable :: cym_str
+  character(len=:), allocatable :: interp_str
   integer :: max_mode, iter, i, m, rst_timestep
   integer, dimension(2,2) :: gc_num
   real :: k0 
@@ -104,10 +105,20 @@ subroutine init_sim_lasers( this, input, opts )
   call input%get( 'simulation.read_restart', read_rst )
   call kwargs%append( 'iter', 0 )
   call kwargs%append( 'k0', 10.0 )
-
-  gc_num(:,1) = (/1, 1/)
-  gc_num(:,2) = (/2, 1/)
-
+
+  call input%get ( 'simulation.interpolation', interp_str)
+  select case ( trim(interp_str) )
+  case ( 'linear' )
+    gc_num(:,1) = (/1, 1/)
+    gc_num(:,2) = (/2, 1/)
+  case ( 'quadratic' )
+    gc_num(:,1) = (/1, 2/)
+    gc_num(:,2) = (/2, 1/)
+  case default
+    call write_err( 'Invalid interpolation type! Only "linear" and "quadratic" are supported &
+      &currently.' )
+  end select
+
   do i = 1, this%num_lasers
     call input%get( 'laser(' // num2str(i) // ').iteration', iter )
     call input%get( 'laser(' // num2str(i) // ').k0', k0 )

source/sim_plasma_class.f03

@@ -88,8 +88,11 @@ subroutine init_sim_plasma( this, input, opts, s )
   select case ( trim(str) )
   case ( 'linear' )
     ps = p_ps_linear
+  case ( 'quadratic' )
+    ps = p_ps_quadratic
+    call write_stdout('Attention! Quadratic interpolation for plasma not finised!')
   case default
-    call write_err( 'Invalid interpolation type! Only "linear" are supported currently.' )
+    call write_err( 'Invalid interpolation type! Only "linear" and "quadratic" are supported currently.' )
   end select
 
   ! initialize profiles of species and neutrals
@@ -170,7 +173,7 @@ subroutine init_sim_plasma( this, input, opts, s )
       call input%get( 'neutrals('//num2str(i)//').smooth_order', sm_ord )
     endif
 
-    call this%neut(i)%new( opts, this%pf_neut(i), max_mode, elem, ion_max, &
+    call this%neut(i)%new( opts, this%pf_neut(i), ps, max_mode, elem, ion_max, &
       qbm, omega_p, s, push_type, sm_ord )
 
   enddo