sll_m_maxwell_1d_trafo.F90

Go to the documentation of this file.

 
 
module sll_m_maxwell_1d_trafo
  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#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_mapping_3d, only: &
       sll_t_mapping_3d
 
  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_multiply_g_1d, &
       sll_s_multiply_gt_1d
 
  use sll_m_spline_fem_utilities_sparse, only : &
       sll_s_spline_fem_mass1d_full, &
       sll_s_spline_fem_mixedmass1d_full
 
  implicit none
 
  public :: &
       sll_t_maxwell_1d_trafo
 
  private
  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
  type, extends(sll_c_maxwell_1d_base) :: sll_t_maxwell_1d_trafo
 
     sll_real64, allocatable :: work1(:)   
     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
     type( sll_t_linear_operator_schur_eb_1d ) :: linear_op_schur_eb
     type( sll_t_linear_solver_mgmres )        :: linear_solver_schur_eb
     type(sll_t_mapping_3d), pointer :: map
     sll_real64 :: force_sign = 1._f64 
 
   contains
 
     procedure :: &
          compute_e_from_b => sll_s_compute_e_from_b_1d_trafo 
     procedure :: &
          compute_b_from_e => sll_s_compute_b_from_e_1d_trafo 
     procedure :: &
          compute_curl_part => sll_s_compute_curl_part_1d_trafo 
     procedure :: &
          compute_e_from_rho => sll_s_compute_e_from_rho_1d_trafo 
     procedure :: &
          compute_rho_from_e => compute_rho_from_e_1d_trafo 
     procedure :: &
          compute_e_from_j => compute_e_from_j_1d_trafo 
     procedure :: &
          compute_phi_from_rho => compute_phi_from_rho_1d_trafo 
     procedure :: &
          compute_phi_from_j => compute_phi_from_j_1d_trafo 
 
     procedure :: &
          compute_rhs_from_function => sll_s_compute_rhs_trafo 
     procedure :: &
          l2projection => l2projection_1d_trafo 
     procedure :: &
          l2norm_squared => l2norm_squared_1d_trafo 
     procedure :: &
          inner_product => inner_product_1d_trafo 
     procedure :: &
          init => init_1d_trafo 
     procedure :: &
          init_from_file => init_from_file_1d_trafo 
     procedure :: &
          free => free_1d_trafo 
     procedure :: &
          multiply_g
     procedure :: &
          multiply_gt
     procedure :: &
          multiply_mass => multiply_mass_1d_trafo 
     procedure :: &
          invert_mass => invert_mass_1d_trafo 
 
  end type sll_t_maxwell_1d_trafo
 
contains
 
  subroutine sll_s_compute_e_from_b_1d_trafo(self, delta_t, field_in, field_out)
    class(sll_t_maxwell_1d_trafo) :: 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%work1 )
    self%work2=0._f64
    call self%linear_solver_mass0%solve ( self%work1, self%work2 )
    ! Update ey from self value
    field_out =  field_out + delta_t*self%work2
 
  end subroutine sll_s_compute_e_from_b_1d_trafo
 
 
  subroutine sll_s_compute_b_from_e_1d_trafo(self, delta_t, field_in, field_out)
    class(sll_t_maxwell_1d_trafo) :: 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%work1)
    ! Update Bz from self value
    field_out = field_out - delta_t * self%work1
 
  end subroutine sll_s_compute_b_from_e_1d_trafo
 
 
  subroutine sll_s_compute_curl_part_1d_trafo( self, delta_t, efield, bfield, betar )
    class(sll_t_maxwell_1d_trafo) :: 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%work1 = 0._f64
    self%work2 = 0._f64
 
    ! Compute D^T M2 b
    call self%multiply_mass( bfield, self%work1, self%s_deg_1 )
    call self%multiply_gt( self%work1, self%work2 )
 
    self%linear_op_schur_eb%sign = -delta_t**2*0.25_f64*factor
    call self%linear_op_schur_eb%dot( efield, self%work1 )
    self%work1 = self%work1 + 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%work1, 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_trafo
 
 
  subroutine sll_s_compute_e_from_rho_1d_trafo(self, field_in, field_out )       
    class(sll_t_maxwell_1d_trafo) :: self
    sll_real64, intent( in    )   :: field_in(:)  
    sll_real64, intent(   out )   :: field_out(:) 
 
    call self%poisson_solver%solve( field_in, self%work1 )
    call multiply_g( self, self%work1, field_out )
    field_out = -field_out
 
  end subroutine sll_s_compute_e_from_rho_1d_trafo
 
 
  subroutine compute_rho_from_e_1d_trafo(self, field_in, field_out)
    class(sll_t_maxwell_1d_trafo) :: self
    sll_real64, intent( in    )   :: field_in(:)  
    sll_real64, intent(   out )   :: field_out(:) 
 
    call self%multiply_mass( field_in, self%work1, self%s_deg_1 )
    call multiply_gt( self, self%work1, field_out )
    field_out = - self%force_sign * field_out
 
  end subroutine compute_rho_from_e_1d_trafo
 
 
  subroutine compute_e_from_j_1d_trafo(self, current, component, E)
    class(sll_t_maxwell_1d_trafo)         :: 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%work1  )
 
    elseif (component == 2) then
       call self%linear_solver_mass0%solve( current, self%work1  )
    else
       print*, 'Component ', component, 'not implemented in compute_E_from_j_1d_trafo.'
    end if
 
    ! Update the electric field and scale
    e = e - self%force_sign * self%work1
 
  end subroutine compute_e_from_j_1d_trafo
 
 
  subroutine compute_phi_from_rho_1d_trafo( self, in, phi, efield )
    class(sll_t_maxwell_1d_trafo)              :: 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_trafo
 
 
  subroutine compute_phi_from_j_1d_trafo( self, in, phi, efield )
    class(sll_t_maxwell_1d_trafo)              :: 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%work1 )
 
    ! Compute phi by inverting the mass matrix
    call self%invert_mass( self%work1, self%work2, self%s_deg_0 )
    phi = phi + self%work2
 
    ! 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_trafo
 
 
  subroutine sll_s_compute_rhs_trafo(self, func, degree, coefs_dofs)
    class(sll_t_maxwell_1d_trafo)       :: 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, c
    sll_real64, dimension(2,degree+1) :: xw_gauss
    sll_real64, dimension(degree+1,degree+1) :: bspl
    sll_real64 :: jacobian
 
    ! 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))
    end do
 
    if( degree == self%s_deg_0 ) then
       ! 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
                c = self%delta_x * (xw_gauss(1,k) + real(i + j - 2 - ((i + j - 2)/self%n_cells) * self%n_cells,f64))
                jacobian= self%map%jacobian( [[c, 0.0_f64, 0._f64]])
                coef = coef + xw_gauss(2,k)*func(self%map%get_x1( [c, 0._f64, 0._f64])) * bspl(degree+2-j,k)*abs(jacobian)
             enddo
          enddo
          ! rescale by cell size
          coefs_dofs(i) = coef*self%delta_x
       enddo
    else
       ! 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
                c = self%delta_x * (xw_gauss(1,k) + real(i + j - 2 - ((i + j - 2)/self%n_cells) * self%n_cells,f64))
                coef = coef + xw_gauss(2,k)*func(self%map%get_x1( [c, 0._f64, 0._f64])) * bspl(degree+2-j,k)* sign( 1._f64, self%map%jacobian( [c, 0._f64, 0._f64] ) )
             enddo
          enddo
          ! rescale by cell size
          coefs_dofs(i) = coef*self%delta_x
       enddo
    endif
 
  end subroutine sll_s_compute_rhs_trafo
 
 
  subroutine l2projection_1d_trafo(self, func, degree, coefs_dofs)
    class(sll_t_maxwell_1d_trafo) :: self
    procedure(sll_i_function_1d_real64) :: func
    sll_int32,  intent( in  ) :: degree 
    sll_real64, intent( out ) :: coefs_dofs(:)  
 
    ! Compute right-hand-side
    call self%compute_rhs_from_function( func, degree, self%work1)
 
    ! Multiply by inverse mass matrix 
    if (degree == self%s_deg_0) then
       call self%linear_solver_mass0%solve ( self%work1, coefs_dofs )
 
    elseif  (degree == self%s_deg_1) then
       call self%linear_solver_mass1%solve ( self%work1, coefs_dofs )
    else
       print*, 'degree ', degree, 'not availlable in maxwell_1d_trafo object' 
    endif
  end subroutine l2projection_1d_trafo
 
 
  function l2norm_squared_1d_trafo(self, coefs_dofs, degree) result (r)
    class(sll_t_maxwell_1d_trafo) :: self
    sll_real64 :: coefs_dofs(:) 
    sll_int32  :: degree 
    sll_real64 :: r 
 
    ! Multiply coefficients by mass matrix (use diagonalization FFT and mass matrix eigenvalues)
    if (degree == self%s_deg_0 ) then
 
       call self%multiply_mass( coefs_dofs, self%work1, self%s_deg_0 )
 
    elseif (degree == self%s_deg_1) then
 
       call self%multiply_mass( coefs_dofs, self%work1, self%s_deg_1 )
 
    end if
    ! Multiply by the coefficients from the left (inner product)
    r = sum(coefs_dofs*self%work1)
 
  end function l2norm_squared_1d_trafo
 
 
  function inner_product_1d_trafo(self, coefs1_dofs, coefs2_dofs, degree, degree2) result (r)
    class(sll_t_maxwell_1d_trafo) :: self
    sll_real64 :: coefs1_dofs(:) 
    sll_real64 :: coefs2_dofs(:) 
    sll_int32  :: degree 
    sll_int32, optional  :: degree2 
    sll_real64 :: r 
 
    ! Multiply coefficients by mass matrix 
    if (present(degree2)) then
       if (degree == degree2) then
          if (degree == self%s_deg_0 ) then
 
             call self%multiply_mass( coefs2_dofs, self%work1, self%s_deg_0 )
 
          elseif (degree == self%s_deg_1) then
 
             call self%multiply_mass( coefs2_dofs, self%work1, self%s_deg_1 )
 
          end if
          ! Multiply by the coefficients from the left (inner product)
          r = sum(coefs1_dofs*self%work1)
       else
          if (degree == self%s_deg_0) then
             call self%multiply_mass( coefs2_dofs, self%work1, 10 )
             ! Multiply by the coefficients from the left (inner product)
             r = sum(coefs1_dofs*self%work1)
          else
             call self%multiply_mass( coefs1_dofs, self%work1, 10 )
             ! Multiply by the coefficients from the left (inner product)
             r = sum(coefs2_dofs*self%work1)
          end if
       end if
    else
       if (degree == self%s_deg_0 ) then
 
          call self%multiply_mass( coefs2_dofs, self%work1, self%s_deg_0 )
 
       elseif (degree == self%s_deg_1) then
 
          call self%multiply_mass( coefs2_dofs, self%work1, self%s_deg_1 )
 
       end if
       ! Multiply by the coefficients from the left (inner product)
       r = sum(coefs1_dofs*self%work1)
    end if
  end function inner_product_1d_trafo
 
 
  subroutine init_1d_trafo( self, domain, n_dofs, s_deg_0, map, mass_tolerance, poisson_tolerance, solver_tolerance, force_sign )
    class(sll_t_maxwell_1d_trafo), intent(out) :: self
    sll_real64, intent(in) :: domain(2)     
    sll_int32, intent(in) :: n_dofs  
    sll_int32, intent(in) :: s_deg_0 
    type(sll_t_mapping_3d), target, intent(inout) :: map
    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 
    ! local variables
    sll_real64, allocatable :: nullspace(:,:)
    sll_int32 :: ierr
 
    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
 
    self%n_dofs0 = n_dofs
    self%n_dofs1 = n_dofs
    self%n_cells = n_dofs
    self%Lx =  map%params(1)
    self%delta_x = 1._f64 / real(n_dofs,f64)
    self%s_deg_0 = s_deg_0
    self%s_deg_1 = s_deg_0 - 1
    self%map=>map
 
    sll_allocate(self%work1(n_dofs),ierr)
    sll_allocate(self%work2(n_dofs),ierr)
 
    ! Sparse matrices
    call sll_s_spline_fem_mass1d_full ( n_dofs, self%s_deg_0, self%mass0, profile_0 )
    call sll_s_spline_fem_mass1d_full ( n_dofs, self%s_deg_1, self%mass1, profile_1 )
 
    call self%linear_solver_mass0%create( self%mass0)
    self%linear_solver_mass0%atol = self%mass_solver_tolerance
 
    call self%linear_solver_mass1%create( self%mass1)
    self%linear_solver_mass1%atol = self%mass_solver_tolerance /self%Lx
    !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_mixedmass1d_full( n_dofs, self%s_deg_0, self%s_deg_1, self%mixed_mass, profile_m0 )
 
    ! 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.
 
  contains
    function profile_0(x)
      sll_real64 :: profile_0
      sll_real64, intent(in) :: x
 
      profile_0 = map%jacobian( [x, 0._f64, 0._f64] )
 
    end function profile_0
 
    function profile_1(x)
      sll_real64 :: profile_1
      sll_real64, intent(in) :: x
 
      profile_1 = 1._f64/map%jacobian( [x, 0._f64, 0._f64] )
 
    end function profile_1
 
    function profile_m0(x)
      sll_real64 :: profile_m0
      sll_real64, intent(in) :: x
 
      profile_m0 = 1._f64
    end function profile_m0
 
  end subroutine init_1d_trafo
 
 
  subroutine init_from_file_1d_trafo( self, domain, n_dofs, s_deg_0, map, nml_file, force_sign )
    class(sll_t_maxwell_1d_trafo), intent(out) :: self
    sll_real64, intent(in) :: domain(2)     
    sll_int32, intent(in) :: n_dofs  
    sll_int32, intent(in) :: s_deg_0 
    type(sll_t_mapping_3d), target, intent(inout) :: map
    character(len=*), intent(in) :: nml_file
    sll_real64, intent(in), optional :: force_sign 
    ! local variables
    sll_int32 :: input_file
    sll_int32 :: io_stat, io_stat0, io_stat1, rank, file_id
    character(len=256) :: file_prefix
    sll_real64 :: mass_tolerance
    sll_real64 :: poisson_tolerance
    sll_real64 :: maxwell_tolerance
 
    namelist /output/ file_prefix
    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
 
    ! 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_trafo: 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, map )
    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 .and. io_stat1 /= 0) then
          if (rank == 0 ) then
             print*, 'sll_m_maxwell_1d_trafo: Input parameter 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, map )
       else if (io_stat /= 0 .and. io_stat1 == 0) then
          call self%init( domain, n_dofs, s_deg_0, map, solver_tolerance=maxwell_tolerance )
       else if (io_stat == 0 .and. io_stat1 /= 0) then
          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, map, mass_tolerance, poisson_tolerance )
       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, map, mass_tolerance, poisson_tolerance, solver_tolerance=maxwell_tolerance )
       end if
       close(input_file)
    end if
 
 
  end subroutine init_from_file_1d_trafo
 
 
  subroutine free_1d_trafo(self)
    class(sll_t_maxwell_1d_trafo) :: self
 
    deallocate(self%work1)
    deallocate(self%work2)
    self%map=> null()
    call self%linear_solver_mass0%free()
    call self%linear_solver_mass1%free()
    call self%poisson_solver%free()
    call self%poisson_operator%free()
    call self%poisson_matrix%free()
    call self%linear_solver_schur_eb%free()
    call self%linear_op_schur_eb%free()
 
  end subroutine free_1d_trafo
 
 
  subroutine multiply_g( self,  in, out)
    class(sll_t_maxwell_1d_trafo), 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_trafo), 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_trafo( self,  in, out, degree)
    class(sll_t_maxwell_1d_trafo), intent(inout) :: self
    sll_real64, intent(in)  :: in(:) 
    sll_real64, intent(out) :: out(:) 
    sll_int32, intent(in)  :: degree 
 
    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_trafo: multiply mass for other form not yet implemented'
       stop
    end if
 
  end subroutine multiply_mass_1d_trafo
 
 
  subroutine invert_mass_1d_trafo( self,  in, out, degree)
    class(sll_t_maxwell_1d_trafo), intent(inout) :: self
    sll_real64, intent(in)  :: in(:) 
    sll_real64, intent(out) :: out(:) 
    sll_int32, intent(in)  :: degree 
 
    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_trafo
 
 
end module sll_m_maxwell_1d_trafo