sll_m_maxwell_1d_fem_sm.F90

Go to the documentation of this file.

 
 
module sll_m_maxwell_1d_fem_sm
  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#include "sll_assert.h"
#include "sll_memory.h"
#include "sll_working_precision.h"
 
  use sll_m_collective, only: &
       sll_f_get_collective_rank, &
       sll_v_world_collective
 
  use sll_m_gauss_legendre_integration, only: &
       sll_f_gauss_legendre_points_and_weights
 
  use sll_m_linear_operator_penalized, only : &
       sll_t_linear_operator_penalized
 
  use sll_m_linear_operator_poisson_1d, only : &
       sll_t_linear_operator_poisson_1d
 
  use sll_m_linear_operator_schur_eb_1d, only : &
       sll_t_linear_operator_schur_eb_1d
 
  use sll_m_linear_solver_cg, only : &
       sll_t_linear_solver_cg
 
  use sll_m_linear_solver_mgmres, only : &
       sll_t_linear_solver_mgmres
 
  use sll_m_low_level_bsplines, only: &
       sll_s_uniform_bsplines_eval_basis
 
  use sll_m_matrix_csr, only: &
       sll_t_matrix_csr
 
  use sll_m_maxwell_1d_base, only: &
       sll_i_function_1d_real64, &
       sll_c_maxwell_1d_base
 
  use sll_m_spline_fem_utilities, only : &
       sll_s_spline_fem_mass_line, &
       sll_s_spline_fem_mixedmass_line, &
       sll_s_multiply_g_1d, &
       sll_s_multiply_gt_1d
 
  use sll_m_spline_fem_utilities_sparse, only : &
       sll_s_spline_fem_mass1d, &
       sll_s_spline_fem_mixedmass1d
 
 
  implicit none
 
  public :: &
       sll_t_maxwell_1d_fem_sm
 
  private
  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  type, extends(sll_c_maxwell_1d_base) :: sll_t_maxwell_1d_fem_sm
 
     type( sll_t_linear_operator_schur_eb_1d ) :: linear_op_schur_eb
     type( sll_t_linear_solver_mgmres )        :: linear_solver_schur_eb
 
     sll_real64, allocatable :: wsave(:) 
     sll_real64, allocatable :: work(:)  
     sll_real64, allocatable :: work2(:)  
     type(sll_t_matrix_csr)  :: mass0
     type(sll_t_matrix_csr)  :: mass1
     type(sll_t_matrix_csr)  :: mixed_mass
     type(sll_t_linear_solver_cg) :: linear_solver_mass0
     type(sll_t_linear_solver_cg) :: linear_solver_mass1
 
     type(sll_t_linear_operator_poisson_1d) :: poisson_matrix
     type(sll_t_linear_operator_penalized) :: poisson_operator
     type(sll_t_linear_solver_cg) :: poisson_solver
     sll_real64 :: force_sign = 1._f64 
     logical ::  adiabatic_electrons
 
   contains
     procedure :: &
          compute_e_from_b => sll_s_compute_e_from_b_1d_fem_sm 
     procedure :: &
          compute_b_from_e => sll_s_compute_b_from_e_1d_fem_sm 
     procedure :: &
          compute_curl_part => sll_s_compute_curl_part_1d_fem_sm 
     procedure :: &
          compute_e_from_rho => sll_s_compute_e_from_rho_1d_fem_sm 
     procedure :: &
          compute_rho_from_e => compute_rho_from_e_1d_fem_sm 
     procedure :: &
          compute_e_from_j => compute_e_from_j_1d_fem_sm 
     procedure :: &
          compute_phi_from_rho => compute_phi_from_rho_1d_fem_sm 
     procedure :: &
          compute_phi_from_j => compute_phi_from_j_1d_fem_sm 
     procedure :: &
          compute_rhs_from_function => sll_s_compute_rhs_fem_sm 
     procedure :: &
          l2projection => l2projection_1d_fem_sm 
     procedure :: &
          l2norm_squared => l2norm_squared_1d_fem_sm 
     procedure :: &
          inner_product => inner_product_1d_fem_sm 
     procedure :: &
          init => init_1d_fem_sm 
     procedure :: &
          init_from_file => init_from_file_1d_fem_sm 
     procedure :: &
          free => free_1d_fem_sm 
     procedure :: &
          multiply_g
     procedure :: &
          multiply_gt
     procedure :: &
          multiply_mass => multiply_mass_1d_fem_sm 
     procedure :: &
          invert_mass => invert_mass_1d_fem_sm 
 
 
  end type sll_t_maxwell_1d_fem_sm
 
contains
 
 
  subroutine sll_s_compute_e_from_b_1d_fem_sm(self, delta_t, field_in, field_out)
    class(sll_t_maxwell_1d_fem_sm) :: self
    sll_real64, intent(in)     :: delta_t   
    sll_real64, intent(in)     :: field_in(:)  
    sll_real64, intent(inout)  :: field_out(:)  
 
    call self%multiply_mass( field_in, self%work2, self%s_deg_1 )
    call self%multiply_gt( self%work2, self%work )
    self%work2=0._f64
    !call self%linear_solver_mass0%set_guess( field_in )
    call self%linear_solver_mass0%solve ( self%work, self%work2 )
    ! Update bz from self value
    field_out =  field_out + delta_t*self%work2
  end subroutine sll_s_compute_e_from_b_1d_fem_sm
 
 
  subroutine sll_s_compute_b_from_e_1d_fem_sm(self, delta_t, field_in, field_out)
    class(sll_t_maxwell_1d_fem_sm)  :: self
    sll_real64, intent(in)     :: delta_t 
    sll_real64, intent(in)     :: field_in(:)  ! ey
    sll_real64, intent(inout)  :: field_out(:) ! bz
 
    call self%multiply_g(field_in, self%work)
    ! Update Bz from self value
    field_out = field_out - delta_t * self%work
 
  end subroutine sll_s_compute_b_from_e_1d_fem_sm
 
 
  subroutine sll_s_compute_curl_part_1d_fem_sm( self, delta_t, efield, bfield, betar )
    class(sll_t_maxwell_1d_fem_sm) :: self
    sll_real64, intent(in)     :: delta_t   
    sll_real64, intent(inout)  :: efield(:)  
    sll_real64, intent(inout)  :: bfield(:)  
    sll_real64, optional       :: betar      
    !local variables
    sll_real64 :: factor
 
    if( present(betar) ) then
       factor = betar
    else
       factor = 1._f64
    end if
 
    self%work = 0._f64
    self%work2 = 0._f64
 
    ! Compute D^T M2 b
    call self%multiply_mass( bfield, self%work, self%s_deg_1 )
    call self%multiply_gt( self%work, self%work2 )
 
    self%linear_op_schur_eb%sign = -delta_t**2*0.25_f64*factor
    call self%linear_op_schur_eb%dot( efield, self%work )
    self%work = self%work + delta_t*factor*self%work2
 
    ! Save efield dofs from previous time step for B field update
    self%work2 = efield
 
    ! Invert Schur complement matrix
    self%linear_op_schur_eb%sign = delta_t**2*0.25_f64*factor
    call self%linear_solver_schur_eb%set_guess( efield )
    call self%linear_solver_schur_eb%solve( self%work, efield )
 
    ! Update B field
    self%work2 = self%work2 + efield
    call self%compute_B_from_E( delta_t*0.5_f64, self%work2, bfield )
 
 
  end subroutine sll_s_compute_curl_part_1d_fem_sm
 
 
  subroutine sll_s_compute_e_from_rho_1d_fem_sm(self, field_in, field_out )       
    class(sll_t_maxwell_1d_fem_sm) :: self
    sll_real64, intent(in)  :: field_in(:) 
    sll_real64, intent(out) :: field_out(:) 
 
    self%wsave = self%force_sign * field_in
    call self%poisson_solver%solve( self%wsave, self%work )
    call multiply_g( self, self%work, field_out )
    field_out = -field_out
 
  end subroutine sll_s_compute_e_from_rho_1d_fem_sm
 
 
  subroutine compute_rho_from_e_1d_fem_sm(self, field_in, field_out)
    class(sll_t_maxwell_1d_fem_sm)  :: self
    sll_real64, intent(in)     :: field_in(:)  
    sll_real64, intent(out)  :: field_out(:) 
 
    call self%multiply_mass( field_in, self%work, self%s_deg_1 )
    call multiply_gt( self, self%work, field_out )
    field_out = - self%force_sign * field_out
 
  end subroutine compute_rho_from_e_1d_fem_sm
 
 
  subroutine compute_e_from_j_1d_fem_sm(self, current, component, E)
    class(sll_t_maxwell_1d_fem_sm)             :: self
    sll_real64,dimension(:),intent(in)    :: current 
    sll_int32, intent(in)                 :: component 
    sll_real64,dimension(:),intent(inout) :: e 
 
    ! Multiply by inverse mass matrix  using the eigenvalues of the circulant inverse matrix
    if (component == 1) then        
       call self%linear_solver_mass1%solve( current, self%work  )
 
    elseif (component == 2) then
       call self%linear_solver_mass0%solve( current, self%work  )
    else
       print*, 'Component ', component, 'not implemented in compute_E_from_j_1d_fem_sm.'
    end if
 
    ! Update the electric field and scale
    e = e - self%force_sign * self%work
  end subroutine compute_e_from_j_1d_fem_sm
 
 
  subroutine compute_phi_from_rho_1d_fem_sm( self, in, phi, efield )
    class(sll_t_maxwell_1d_fem_sm)             :: self
    sll_real64, intent(in)                     :: in(:)     
    sll_real64, intent(out)                    :: phi(:)    
    sll_real64, intent(out)                    :: efield(:) 
 
    ! Compute phi by inverting the mass matrix
    call self%invert_mass( in, phi, self%s_deg_0 )
 
    ! Compute the degrees of freedom of the electric field as -G phi
    call multiply_g( self, phi, efield )
    efield = -efield
 
  end subroutine compute_phi_from_rho_1d_fem_sm
 
 
  subroutine compute_phi_from_j_1d_fem_sm( self, in, phi, efield )
    class(sll_t_maxwell_1d_fem_sm)             :: self
    sll_real64, intent(in)                     :: in(:)     
    sll_real64, intent(out)                    :: phi(:)    
    sll_real64, intent(out)                    :: efield(:) 
 
    ! Compute divergence of the current (G^T current) (assuming a 1v model)
    call multiply_gt( self, in, self%work )
 
    ! Compute phi by inverting the mass matrix
    call self%invert_mass( self%work, self%wsave, self%s_deg_0 )
    phi = phi + self%wsave
 
    ! Compute the degrees of freedom of the electric field as -G phi
    call multiply_g( self, phi, efield )
    efield = -efield
 
  end subroutine compute_phi_from_j_1d_fem_sm
 
 
  subroutine sll_s_compute_rhs_fem_sm(self, func, degree, coefs_dofs)
    class(sll_t_maxwell_1d_fem_sm)      :: self
    procedure(sll_i_function_1d_real64) :: func
    sll_int32, intent(in) :: degree 
    sll_real64, intent(out) :: coefs_dofs(:)  
    ! local variables
    sll_int32 :: i,j,k
    sll_real64 :: coef
    sll_real64, dimension(2,degree+1) :: xw_gauss
    sll_real64, dimension(degree+1,degree+1) :: bspl
 
    ! take enough Gauss points so that projection is exact for splines of degree deg
    ! rescale on [0,1] for compatibility with B-splines
    xw_gauss = sll_f_gauss_legendre_points_and_weights(degree+1, 0.0_f64, 1.0_f64)
    ! Compute bsplines at gauss_points
    do k=1,degree+1
       call sll_s_uniform_bsplines_eval_basis(degree,xw_gauss(1,k), bspl(:,k))
       !print*, 'bs', bspl(k,:)
    end do
 
    ! Compute coefs_dofs = int f(x)N_i(x) 
    do i = 1, self%n_cells
       coef=0.0_f64
       ! loop over support of B spline
       do j = 1, degree+1
          ! loop over Gauss points
          do k=1, degree+1
             coef = coef + xw_gauss(2,k)*func(self%delta_x * (xw_gauss(1,k) + real(i + j - 2 - ((i + j - 2)/self%n_cells) * self%n_cells,f64))) * bspl(degree+2-j,k)
          end do
       end do
       ! rescale by cell size
       coefs_dofs(i) = coef*self%delta_x
    end do
 
  end subroutine sll_s_compute_rhs_fem_sm
 
 
  subroutine l2projection_1d_fem_sm(self, func, degree, coefs_dofs)
    class(sll_t_maxwell_1d_fem_sm) :: self
    procedure(sll_i_function_1d_real64) :: func
    sll_int32, intent(in) :: degree 
    sll_real64, intent(out) :: coefs_dofs(:)  
 
    ! Compute right-hand-side
    call sll_s_compute_rhs_fem_sm(self, func, degree, self%work)
 
    ! Multiply by inverse mass matrix (! complex numbers stored in real array with fftpack ordering)
    if (degree == self%s_deg_0) then
       call self%linear_solver_mass0%solve ( self%work, coefs_dofs )
 
    elseif  (degree == self%s_deg_0-1) then
       call self%linear_solver_mass1%solve ( self%work, coefs_dofs )
    else
       print*, 'degree ', degree, 'not availlable in maxwell_1d_fem_sm object' 
    endif
 
  end subroutine l2projection_1d_fem_sm
 
 
  function l2norm_squared_1d_fem_sm(self, coefs_dofs, degree) result (r)
    class(sll_t_maxwell_1d_fem_sm) :: self
    sll_real64 :: coefs_dofs(:) 
    sll_int32  :: degree 
    sll_real64 :: r 
 
    ! Multiply coefficients by mass matrix 
    if (degree == self%s_deg_0 ) then
       call self%multiply_mass( coefs_dofs, self%work, self%s_deg_0 )
 
    elseif (degree == self%s_deg_1) then
       call self%multiply_mass( coefs_dofs, self%work, self%s_deg_1 )
 
    end if
    ! Multiply by the coefficients from the left (inner product)
    r = sum(coefs_dofs*self%work)
  end function l2norm_squared_1d_fem_sm
 
 
  function inner_product_1d_fem_sm(self, coefs1_dofs, coefs2_dofs, degree, degree2) result (r)
    class(sll_t_maxwell_1d_fem_sm) :: self
    sll_real64 :: coefs1_dofs(:) 
    sll_real64 :: coefs2_dofs(:) 
    sll_int32  :: degree 
    sll_int32, optional  :: degree2 
    sll_real64 :: r 
 
    if (present(degree2)) then
       if (degree == degree2) then
          if (degree == self%s_deg_0 ) then
 
             !call self%mass0%dot ( coefs2_dofs, self%work )
             call self%multiply_mass( coefs2_dofs, self%work, self%s_deg_0 )
 
          elseif (degree == self%s_deg_1) then
 
             call self%multiply_mass( coefs2_dofs, self%work, self%s_deg_1 )
             !call self%mass1%dot ( coefs2_dofs, self%work )
 
          end if
 
          ! Multiply by the coefficients from the left (inner product)
          r = sum(coefs1_dofs*self%work)
       else
          if (degree == self%s_deg_0) then
             call self%multiply_mass( coefs2_dofs, self%work, 10 )
             ! Multiply by the coefficients from the left (inner product)
             r = sum(coefs1_dofs*self%work)
          else
             call self%multiply_mass( coefs1_dofs, self%work, 10 )
             ! Multiply by the coefficients from the left (inner product)
             r = sum(coefs2_dofs*self%work)
          end if
       end if
    else
       if (degree == self%s_deg_0 ) then
          call self%multiply_mass( coefs2_dofs, self%work, self%s_deg_0 )
 
       elseif (degree == self%s_deg_1) then
          call self%multiply_mass( coefs2_dofs, self%work, self%s_deg_1 )
 
       end if
       ! Multiply by the coefficients from the left (inner product)
       r = sum(coefs1_dofs*self%work)
    end if
  end function inner_product_1d_fem_sm
 
 
  subroutine init_1d_fem_sm( self, domain, n_dofs, s_deg_0, delta_t,  mass_tolerance, poisson_tolerance, solver_tolerance, force_sign, adiabatic_electrons )
    class(sll_t_maxwell_1d_fem_sm), intent(out) :: self
    sll_real64, intent(in) :: domain(2)     
    sll_int32, intent(in) :: n_dofs  
    !sll_real64 :: delta_x ! cell size
    sll_int32, intent(in) :: s_deg_0 
    sll_real64, intent(in) :: delta_t 
    sll_real64, intent(in), optional :: mass_tolerance 
    sll_real64, intent(in), optional :: poisson_tolerance 
    sll_real64, intent(in), optional :: solver_tolerance 
    sll_real64, intent(in), optional :: force_sign 
    logical, intent(in), optional :: adiabatic_electrons
    ! local variables
    sll_int32 :: ierr
    sll_real64 :: mass_line_0(s_deg_0+1), mass_line_1(s_deg_0), mass_line_mixed(s_deg_0*2)
    sll_real64, allocatable :: nullspace(:,:)
    sll_real64 :: t_e = 10000._f64
 
    if (present( mass_tolerance) ) then
       self%mass_solver_tolerance = mass_tolerance
    else
       self%mass_solver_tolerance = 1d-12
    end if
    if (present( poisson_tolerance) ) then
       self%poisson_solver_tolerance = poisson_tolerance
    else
       self%poisson_solver_tolerance = 1d-12
    end if
 
    if (present( solver_tolerance) ) then
       self%solver_tolerance = solver_tolerance
    else
       self%solver_tolerance = 1d-12
    end if
 
    if (present( force_sign) ) then
       self%force_sign = force_sign
    end if
 
    if (present( adiabatic_electrons) ) then
       self%adiabatic_electrons = adiabatic_electrons
    end if
 
    self%n_dofs0 = n_dofs
    self%n_dofs1 = n_dofs
    self%n_cells = n_dofs
    self%Lx = domain(2) - domain(1)
    self%delta_x = self%Lx /real(n_dofs,f64)
    self%s_deg_0 = s_deg_0
    self%s_deg_1 = s_deg_0 - 1
 
 
    sll_allocate(self%work(n_dofs),ierr)
    sll_allocate(self%work2(n_dofs),ierr)
    sll_allocate(self%wsave(n_dofs),ierr)
 
    ! Sparse matrices
    call sll_s_spline_fem_mass_line ( self%s_deg_0, mass_line_0 )
    call sll_s_spline_fem_mass_line ( self%s_deg_1, mass_line_1 )
    !scale with delta_x=L/n_dofs
    if(self%adiabatic_electrons) then
       mass_line_0= mass_line_0*self%delta_x/t_e
    else
       mass_line_0= mass_line_0*self%delta_x
    end if
    mass_line_1= mass_line_1*self%delta_x
    call sll_s_spline_fem_mass1d( n_dofs, self%s_deg_0, mass_line_0, self%mass0 )
    call sll_s_spline_fem_mass1d( n_dofs, self%s_deg_1, mass_line_1, self%mass1 )
    call self%linear_solver_mass0%create( self%mass0)
    call self%linear_solver_mass1%create( self%mass1)
    self%linear_solver_mass0%atol = self%mass_solver_tolerance
    self%linear_solver_mass1%atol = self%mass_solver_tolerance
    !self%linear_solver_mass0%verbose = .true.
    !self%linear_solver_mass1%verbose = .true.
 
    ! Penalized Poisson operator
    allocate(nullspace(1,1:self%n_cells))
    nullspace(1,:) = 1.0_f64
    call self%poisson_matrix%create( self%mass1, self%n_cells, self%delta_x )
    call self%poisson_operator%create( self%poisson_matrix, vecs=nullspace, &
         n_dim_nullspace=1 )
    call self%poisson_solver%create( self%poisson_operator )
    self%poisson_solver%null_space = .true.
    self%poisson_solver%atol = self%poisson_solver_tolerance
    !self%poisson_solver%verbose = .true.
 
 
    call sll_s_spline_fem_mixedmass_line( self%s_deg_0, mass_line_mixed)
    mass_line_mixed=mass_line_mixed*self%delta_x
    call sll_s_spline_fem_mixedmass1d( n_dofs, self%s_deg_0, mass_line_mixed, self%mixed_mass )
 
    ! Only for Schur complement eb solver
    call self%linear_op_schur_eb%create( self%mass0, self%mass1, self%n_cells, self%delta_x )
    call self%linear_solver_schur_eb%create( self%linear_op_schur_eb )
    self%linear_solver_schur_eb%atol = self%solver_tolerance
    self%linear_solver_schur_eb%rtol = self%solver_tolerance
    !self%linear_solver_schur_eb%verbose = .true.
 
  end subroutine init_1d_fem_sm
 
 
  subroutine init_from_file_1d_fem_sm( self, domain, n_dofs, s_deg_0, nml_file, force_sign, adiabatic_electrons )
    class(sll_t_maxwell_1d_fem_sm), intent(out) :: self
    sll_real64, intent(in) :: domain(2)     
    sll_int32, intent(in) :: n_dofs  
    sll_int32, intent(in) :: s_deg_0 
    character(len=*) :: nml_file
    sll_real64, intent(in), optional :: force_sign 
    logical, intent(in), optional :: adiabatic_electrons
    ! local variables
    sll_int32 :: input_file, rank
    sll_int32 :: io_stat, io_stat0, io_stat1, file_id 
    character(len=256) :: file_prefix
    logical            :: output_fields = .false.
    logical            :: output_particles = .false.
    sll_real64 :: mass_tolerance
    sll_real64 :: poisson_tolerance
    sll_real64 :: maxwell_tolerance
    sll_real64         :: delta_t = 0._f64
    sll_int32          :: n_time_steps
    sll_real64         :: beta
    character(len=256) :: initial_distrib
    character(len=256) :: initial_bfield
    sll_real64         :: charge
    
    namelist /sim_params/ delta_t, n_time_steps, beta, initial_distrib, initial_bfield, charge
    namelist /output/ file_prefix, output_fields, output_particles
    namelist /maxwell_solver/ mass_tolerance, poisson_tolerance
    namelist /time_solver/ maxwell_tolerance 
 
    rank = sll_f_get_collective_rank(sll_v_world_collective)
 
    if( present(force_sign) ) then
       self%force_sign = force_sign
    end if
 
    if (present( adiabatic_electrons) ) then
       self%adiabatic_electrons = adiabatic_electrons
    end if
 
    ! Read in solver tolerance
    open(newunit = input_file, file=nml_file, status='old',iostat=io_stat)
    if (io_stat /= 0) then
       if (rank == 0 ) then
          print*, 'sll_m_maxwell_1d_fem_sm: Input file does not exist. Set default tolerance.'
          open(newunit=file_id, file=trim(file_prefix)//'_parameters_used.dat', position = 'append', status='old', action='write', iostat=io_stat)
          write(file_id, *) 'mass solver tolerance:', 1d-12
          write(file_id, *) 'poisson solver tolerance:', 1d-12
          close(file_id)
       end if
       call self%init( domain, n_dofs, s_deg_0, delta_t*0.5_f64 )
 
    else
       read(input_file, output,iostat=io_stat0)
       read(input_file, maxwell_solver,iostat=io_stat)
       read(input_file, time_solver,iostat=io_stat1)
       if (io_stat /= 0) then
          if (rank == 0 ) then
             print*, 'sll_m_maxwell_1d_fem_sm: Tolerance not given. Set default tolerance.'
             open(newunit=file_id, file=trim(file_prefix)//'_parameters_used.dat', position = 'append', status='old', action='write', iostat=io_stat)
             write(file_id, *) 'mass solver tolerance:', 1d-12
             write(file_id, *) 'poisson solver tolerance:', 1d-12
             close(file_id)
          end if
          call self%init( domain, n_dofs, s_deg_0, delta_t*0.5_f64 )
       else
          if (rank == 0 ) then
             open(newunit=file_id, file=trim(file_prefix)//'_parameters_used.dat', position = 'append', status='old', action='write', iostat=io_stat)
             write(file_id, *) 'mass solver tolerance:', mass_tolerance
             write(file_id, *) 'poisson solver tolerance:', poisson_tolerance
             close(file_id)
          end if
          call self%init( domain, n_dofs, s_deg_0, delta_t*0.5_f64, mass_tolerance, poisson_tolerance, solver_tolerance=maxwell_tolerance )
       end if
       close(input_file)
    end if
 
 
  end subroutine init_from_file_1d_fem_sm
 
 
  subroutine free_1d_fem_sm(self)
    class(sll_t_maxwell_1d_fem_sm) :: self
 
    deallocate(self%wsave)
    deallocate(self%work)
    call self%linear_solver_schur_eb%free()
    call self%linear_op_schur_eb%free()
 
  end subroutine free_1d_fem_sm
 
 
  subroutine multiply_g( self,  in, out)
    class(sll_t_maxwell_1d_fem_sm), intent(in) :: self
    sll_real64, intent(in)  :: in(:)  
    sll_real64, intent(out) :: out(:) 
 
    call sll_s_multiply_g_1d( self%n_cells, self%delta_x, in, out )
 
  end subroutine multiply_g
 
 
  subroutine multiply_gt( self,  in, out)
    class(sll_t_maxwell_1d_fem_sm), intent(in) :: self
    sll_real64, intent(in)  :: in(:)  
    sll_real64, intent(out) :: out(:) 
 
    call sll_s_multiply_gt_1d( self%n_cells, self%delta_x, in, out )
 
  end subroutine multiply_gt
 
 
  subroutine multiply_mass_1d_fem_sm( self,  in, out, degree)
    class(sll_t_maxwell_1d_fem_sm), intent(inout) :: self
    sll_real64, intent(in)  :: in(:) 
    sll_real64, intent(out) :: out(:) 
    sll_int32, intent(in)  :: degree 
 
    ! Multiply coefficients by mass matrix 
    if (degree == self%s_deg_0 ) then
       call self%mass0%dot ( in, out)
    elseif (degree == self%s_deg_1) then
       call self%mass1%dot ( in, out)
    elseif(degree == 10) then
       call self%mixed_mass%dot( in, out )
    else
       print*, 'maxwell_solver_1d_fem_sm: multiply mass for other form not yet implemented'
       stop
    end if
 
  end subroutine multiply_mass_1d_fem_sm
 
 
  subroutine invert_mass_1d_fem_sm( self,  in, out, degree)
    class(sll_t_maxwell_1d_fem_sm), intent(inout) :: self
    sll_real64, intent(in)  :: in(:) 
    sll_real64, intent(out) :: out(:) 
    sll_int32, intent(in)  :: degree 
 
    ! Multiply coefficients by mass matrix 
    if (degree == self%s_deg_0 ) then
 
       call self%linear_solver_mass0%solve( in, out)
 
    elseif (degree == self%s_deg_1) then
 
       call self%linear_solver_mass1%solve ( in, out)
    else
       print*, 'Invert mass for other form not yet implemented'
       stop
    end if
 
  end subroutine invert_mass_1d_fem_sm
 
end module sll_m_maxwell_1d_fem_sm