sll_m_maxwell_clamped_1d_trafo.F90

Go to the documentation of this file.

 
 
module sll_m_maxwell_clamped_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_poisson_clamped_1d, only : &
       sll_t_linear_operator_poisson_clamped_1d
 
  use sll_m_linear_operator_schur_eb_cl_1d, only : &
       sll_t_linear_operator_schur_eb_cl_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_3d_clamped, only: &
       sll_s_multiply_g_clamped_1d, &
       sll_s_multiply_gt_clamped_1d
 
  use sll_m_spline_fem_utilities_sparse, only : &
       sll_s_spline_fem_mass1d_clamped_full, &
        sll_s_spline_fem_mixedmass1d_clamped_full
       
  use sll_m_splines_pp, only: &
       sll_t_spline_pp_1d, &
       sll_s_spline_pp_init_1d, &
       sll_s_spline_pp_free_1d, &
       sll_f_spline_pp_horner_1d
 
  implicit none
 
  public :: &
       sll_t_maxwell_clamped_1d_trafo
 
  private
  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
  type, extends(sll_c_maxwell_1d_base) :: sll_t_maxwell_clamped_1d_trafo
 
     type(sll_t_spline_pp_1d) :: spline0_pp
     type(sll_t_spline_pp_1d) :: spline1_pp
     sll_real64, allocatable :: work1(:)  
     sll_real64, allocatable :: work0(:)  
     sll_real64, allocatable :: work01(:)  
     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_clamped_1d) :: poisson_matrix
     type(sll_t_linear_solver_cg) :: poisson_solver
     type( sll_t_linear_operator_schur_eb_cl_1d ) :: linear_op_schur_eb
     type( sll_t_linear_solver_mgmres )        :: linear_solver_schur_eb
     type(sll_t_mapping_3d), pointer :: map
 
   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_clamped_1d_trafo
 
contains
 
 
  subroutine sll_s_compute_e_from_b_1d_trafo(self, delta_t, field_in, field_out)
    class(sll_t_maxwell_clamped_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%work1, self%s_deg_1 )
    call self%multiply_gt( self%work1, self%work0 )
 
    call self%invert_mass( self%work0, self%work01, self%s_deg_0 )
    ! Update Ey from self value
    field_out = field_out + delta_t*self%work01
 
  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_clamped_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_clamped_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%work0 = 0._f64
    self%work01 = 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%work0 )
 
    self%linear_op_schur_eb%sign = -delta_t**2*0.25_f64*factor
    call self%linear_op_schur_eb%dot( efield, self%work01 )
    self%work01 = self%work01 + delta_t*factor*self%work0
 
    ! Save efield dofs from previous time step for B field update
    self%work0 = 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%work01, efield )
 
    ! Update B field
    self%work0 = self%work0 + efield
    call self%compute_B_from_E( delta_t*0.5_f64, self%work0, 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_clamped_1d_trafo) :: self
    sll_real64, intent( in    )           :: field_in(:)  
    sll_real64, intent(   out )           :: field_out(:) 
 
    call self%poisson_solver%solve( field_in, self%work0 )
    call multiply_g( self, self%work0, 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_clamped_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 = - 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_clamped_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%invert_mass( current, self%work1, self%s_deg_1  )
       e = e -  self%work1
    elseif (component == 2) then
       call self%invert_mass( current, self%work0, self%s_deg_0  )
       e = e -  self%work0
    else
       print*, 'Component ', component, 'not implemented in compute_E_from_j_1d_trafo.'
    end if
 
  end subroutine compute_e_from_j_1d_trafo
 
 
  subroutine compute_phi_from_rho_1d_trafo( self, in, phi, efield )
    class(sll_t_maxwell_clamped_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_clamped_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%work01 )
 
    ! Compute phi by inverting the mass matrix
    call self%invert_mass( self%work01, self%work0, self%s_deg_0 )
    phi = phi + self%work0
 
    ! 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_clamped_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,quad,q
    sll_real64 :: c
    sll_real64, allocatable :: xw_gauss(:,:)
    sll_real64, allocatable :: bspl(:,:)
    sll_real64 :: jacobian
 
    q = 2*degree+1
    allocate( xw_gauss(2,q) )
    allocate( bspl(degree+1,q) )
    ! 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(q, 0.0_f64, 1.0_f64)
 
    if( degree == self%s_deg_0 ) then
       coefs_dofs = 0._f64
       !first degree-1 cells are different
       !loop over cells
       do i = 1, degree-1
          !loop over splines
          do j = 1, degree+1
             ! loop over Gauss points
             do quad = 1, q
                c = self%delta_x * (xw_gauss(1,quad)+ real(i - 1,f64))
                jacobian= self%map%jacobian( [c, 0._f64, 0._f64])
                coefs_dofs(i-1+j) = coefs_dofs(i-1+j) + &
                     self%delta_x* xw_gauss(2,quad) * func(self%map%get_x1( [c, 0._f64, 0._f64])) *&
                     sll_f_spline_pp_horner_1d( degree, self%spline0_pp%poly_coeffs_boundary_left(:,:,i), xw_gauss(1,quad), j)*abs(jacobian)
             enddo
          enddo
       enddo
       bspl=0._f64
       ! Compute bsplines at gauss_points
       do quad=1,q
          do j=1, degree+1
             bspl(j,quad) = sll_f_spline_pp_horner_1d( degree, self%spline0_pp%poly_coeffs, xw_gauss(1,quad), j)
          end do
       end do
 
       ! Compute coefs_dofs = int f(x)N_i(x)
       !loop over cells
       do i = degree, self%n_cells-degree+1
          ! loop over splines
          do j = 1, degree+1
             ! loop over Gauss points
             do quad = 1, q
                c = self%delta_x * (xw_gauss(1,quad)+ real(i - 1,f64))
                jacobian= self%map%jacobian( [c, 0._f64, 0._f64])
                coefs_dofs(i-1+j) = coefs_dofs(i-1+j) + &
                     self%delta_x * xw_gauss(2,quad) * func(self%map%get_x1( [c, 0._f64, 0._f64])) * bspl(j,quad) * abs(jacobian)
             enddo
          enddo
       enddo
 
       !last degree-1 cells are different
       !loop over cells
       do i = self%n_cells-degree+2, self%n_cells
          !loop over splines
          do j = 1, degree+1
             ! loop over Gauss points
             do quad = 1, q
                c = self%delta_x * (xw_gauss(1,quad)+ real(i - 1,f64))
                jacobian= self%map%jacobian( [c, 0._f64, 0._f64])
                coefs_dofs(i-1+j) = coefs_dofs(i-1+j) + &
                     self%delta_x * xw_gauss(2,quad) * func(self%map%get_x1( [c, 0._f64, 0._f64])) *&
                     sll_f_spline_pp_horner_1d( degree, self%spline0_pp%poly_coeffs_boundary_right(:,:,i-self%n_cells+degree-1), xw_gauss(1,quad), j) * abs(jacobian)
             enddo
          enddo
       end do
 
    else if (degree == self%s_deg_1)then
 
       coefs_dofs = 0._f64
       !first degree-1 cells are different
       !loop over cells
       do i = 1, degree-1
          !loop over splines
          do j = 1, degree+1
             ! loop over Gauss points
             do quad = 1, q
                c = self%delta_x * (xw_gauss(1,quad) + real(i - 1,f64))
                coefs_dofs(i-1+j) = coefs_dofs(i-1+j) + &
                     self%delta_x * xw_gauss(2,quad)*func(self%map%get_x1( [c, 0._f64, 0._f64])) *&
                     sll_f_spline_pp_horner_1d( degree, self%spline1_pp%poly_coeffs_boundary_left(:,:,i), xw_gauss(1,quad), j) * sign( 1._f64, self%map%jacobian( [c, 0._f64, 0._f64] ) )
             enddo
          enddo
       enddo
 
       bspl=0._f64
       ! Compute bsplines at gauss_points
       do quad=1, q
          do j=1, degree+1
             bspl(j,quad) = sll_f_spline_pp_horner_1d( degree, self%spline1_pp%poly_coeffs, xw_gauss(1,quad), j)
          end do
       end do
 
       ! Compute coefs_dofs = int f(x)N_i(x)
       !loop over cells
       do i = degree, self%n_cells-degree+1
          !loop over splines
          do j = 1, degree+1
             ! loop over Gauss points
             do quad = 1, q
                c = self%delta_x * (xw_gauss(1,quad)+ real(i - 1,f64))
                coefs_dofs(i-1+j) = coefs_dofs(i-1+j) + &
                     self%delta_x * xw_gauss(2,quad) * func(self%map%get_x1( [c, 0._f64, 0._f64])) * bspl(j,quad) * sign( 1._f64, self%map%jacobian( [c, 0._f64, 0._f64] ) )
             enddo
          enddo
       enddo
 
       !last degree-1 cells are different
       !loop over cells
       do i = self%n_cells-degree+2, self%n_cells
          !loop over splines
          do j = 1, degree+1
             ! loop over Gauss points
             do quad = 1, q
                c = self%delta_x * (xw_gauss(1,quad)+ real(i - 1,f64))
                coefs_dofs(i-1+j) = coefs_dofs(i-1+j) + &
                     self%delta_x * xw_gauss(2,quad) * func(self%map%get_x1( [c, 0._f64, 0._f64])) *&
                     sll_f_spline_pp_horner_1d( degree, self%spline1_pp%poly_coeffs_boundary_right(:,:,i-self%n_cells+degree-1), xw_gauss(1,quad), j) * sign( 1._f64, self%map%jacobian( [c, 0._f64, 0._f64] ) )
             enddo
          enddo
       end do
 
    endif
 
  end subroutine sll_s_compute_rhs_trafo
 
 
  subroutine l2projection_1d_trafo(self, func, degree, coefs_dofs)
    class(sll_t_maxwell_clamped_1d_trafo) :: self
    procedure(sll_i_function_1d_real64) :: func
    sll_int32,  intent( in  ) :: degree 
    sll_real64, intent( out ) :: coefs_dofs(:)  
 
    ! Multiply by inverse mass matrix (! complex numbers stored in real array with fftpack ordering)
    if (degree == self%s_deg_0) then
       ! Compute right-hand-side
       call sll_s_compute_rhs_trafo(self, func, degree, self%work0)
       call self%invert_mass( self%work0, coefs_dofs, self%s_deg_0 )
    elseif  (degree == self%s_deg_1) then
       ! Compute right-hand-side
       call sll_s_compute_rhs_trafo(self, func, degree, self%work1)
       call self%invert_mass( self%work1, coefs_dofs, self%s_deg_1 )
    else
       print*, 'degree ', degree, 'not availlable in maxwell_clamped_1d_trafo object' 
    endif
  end subroutine l2projection_1d_trafo
 
 
  function l2norm_squared_1d_trafo(self, coefs_dofs, degree) result (r)
    class(sll_t_maxwell_clamped_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%work0, self%s_deg_0 )
       ! Multiply by the coefficients from the left (inner product)
       r = sum(coefs_dofs*self%work0)
    elseif (degree == self%s_deg_1) then
       call self%multiply_mass( coefs_dofs, self%work1, self%s_deg_1 )
       ! Multiply by the coefficients from the left (inner product)
       r = sum(coefs_dofs*self%work1)
    end if
 
  end function l2norm_squared_1d_trafo
 
  function inner_product_1d_trafo(self, coefs1_dofs, coefs2_dofs, degree, degree2) result (r)
    class(sll_t_maxwell_clamped_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 (use diagonalization FFT and mass matrix eigenvalues)
    if (present(degree2)) then
       if (degree == degree2) then
          if (degree == self%s_deg_0 ) then
             call self%multiply_mass( coefs2_dofs, self%work0, self%s_deg_0 )
             ! Multiply by the coefficients from the left (inner product)
             r = sum(coefs1_dofs*self%work0)
          elseif (degree == self%s_deg_1) then
             call self%multiply_mass( coefs2_dofs, self%work1, self%s_deg_1 )
             ! Multiply by the coefficients from the left (inner product)
             r = sum(coefs1_dofs*self%work1)
          end if
       else
          if (degree == self%s_deg_0) then
             call self%multiply_mass( coefs2_dofs, self%work01, 10 )
             ! Multiply by the coefficients from the left (inner product)
             r = sum(coefs1_dofs*self%work01)
          else
             call self%multiply_mass( coefs1_dofs, self%work01, 10 )
             ! Multiply by the coefficients from the left (inner product)
             r = sum(coefs2_dofs*self%work01)
          end if
       end if
    else
       if (degree == self%s_deg_0 ) then
          call self%multiply_mass( coefs2_dofs, self%work0, self%s_deg_0 )
          ! Multiply by the coefficients from the left (inner product)
          r = sum(coefs1_dofs*self%work0)
       elseif (degree == self%s_deg_1) then
          call self%multiply_mass( coefs2_dofs, self%work1, self%s_deg_1 )
          ! Multiply by the coefficients from the left (inner product)
          r = sum(coefs1_dofs*self%work1)
       end if
    end if
 
  end function inner_product_1d_trafo
 
 
  subroutine init_1d_trafo( self, domain, n_cells, s_deg_0, boundary, map, mass_tolerance, poisson_tolerance, solver_tolerance )
    class(sll_t_maxwell_clamped_1d_trafo), intent(out) :: self
    sll_real64, intent(in) :: domain(2)     
    sll_int32, intent(in) :: n_cells  
    sll_int32, intent(in) :: s_deg_0 
    sll_int32, intent(in) :: boundary 
    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 
    ! local variables
    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
 
    self%n_cells = n_cells
    self%n_dofs0 = n_cells+s_deg_0
    self%n_dofs1 = n_cells+s_deg_0-1
    self%Lx =  map%params(1)
    self%delta_x = 1._f64 / real(n_cells,f64)
    self%s_deg_0 = s_deg_0
    self%s_deg_1 = s_deg_0 - 1
    self%map=>map
 
    call sll_s_spline_pp_init_1d( self%spline0_pp, self%s_deg_0, n_cells, boundary)
    call sll_s_spline_pp_init_1d( self%spline1_pp, self%s_deg_1, n_cells, boundary)
 
    sll_allocate(self%work1(self%n_dofs1),ierr)
    sll_allocate(self%work0(self%n_dofs0),ierr)
    sll_allocate(self%work01(self%n_dofs0),ierr)
 
    ! Sparse matrices
    call sll_s_spline_fem_mass1d_clamped_full ( n_cells, self%s_deg_0, self%spline0_pp, self%mass0, profile_0 )
    call sll_s_spline_fem_mass1d_clamped_full ( n_cells, self%s_deg_1, self%spline1_pp, 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.
 
    call sll_s_spline_fem_mixedmass1d_clamped_full( n_cells,  self%s_deg_0,  self%s_deg_1, self%mixed_mass, profile_m0, self%spline0_pp, self%spline1_pp )
 
 
    ! Poisson solver
    call self%poisson_matrix%create( self%mass1, self%s_deg_0, self%n_dofs0, self%delta_x )
    call self%poisson_solver%create( self%poisson_matrix )
    self%poisson_solver%atol = self%poisson_solver_tolerance
    !self%poisson_solver%verbose = .true.
    !self%poisson_solver%n_maxiter=40000
 
    ! Only for Schur complement eb solver
    call self%linear_op_schur_eb%create( self%mass0, self%mass1, self%n_dofs0, self%delta_x, self%s_deg_0 )
    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_cells, s_deg_0, boundary, map, nml_file )
    class(sll_t_maxwell_clamped_1d_trafo), intent(out) :: self
    sll_real64, intent(in) :: domain(2)     
    sll_int32, intent(in) :: n_cells  
    sll_int32, intent(in) :: s_deg_0 
    sll_int32, intent(in) :: boundary 
    type(sll_t_mapping_3d), target, intent(inout) :: map
    character(len=*), intent(in) :: nml_file
    ! local variables
    character(len=256) :: file_prefix
    sll_int32 :: input_file
    sll_int32 :: io_stat, io_stat0, io_stat1, rank, file_id
    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)
 
    ! 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_clamped_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_cells, s_deg_0, boundary, 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 .or. io_stat1 /= 0 ) then
          if (rank == 0 ) then
             print*, 'sll_m_maxwell_clamped_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_cells, s_deg_0, boundary, map )
       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_cells, s_deg_0, boundary, map, mass_tolerance, poisson_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_clamped_1d_trafo) :: self
 
    self%map=> null()
    call self%linear_solver_mass0%free()
    call self%linear_solver_mass1%free()
    call self%poisson_solver%free()
    call self%poisson_matrix%free()
    call self%linear_solver_schur_eb%free()
    call self%linear_op_schur_eb%free()
 
    call sll_s_spline_pp_free_1d( self%spline0_pp )
    call sll_s_spline_pp_free_1d( self%spline1_pp )
 
    deallocate(self%work1)
    deallocate(self%work0)
    deallocate(self%work01)
 
  end subroutine free_1d_trafo
 
 
  subroutine multiply_g( self,  in, out)
    class(sll_t_maxwell_clamped_1d_trafo), intent(in) :: self
    sll_real64, intent(in)  :: in(:) 
    sll_real64, intent(out) :: out(:) 
 
    call sll_s_multiply_g_clamped_1d( self%n_dofs0, self%delta_x, self%s_deg_0, in, out)
 
  end subroutine multiply_g
 
 
  subroutine multiply_gt( self,  in, out)
    class(sll_t_maxwell_clamped_1d_trafo), intent(in) :: self
    sll_real64, intent(in)  :: in(:) 
    sll_real64, intent(out) :: out(:) 
 
    call sll_s_multiply_gt_clamped_1d( self%n_dofs0, self%delta_x, self%s_deg_0, in, out)
 
  end subroutine multiply_gt
 
 
  subroutine multiply_mass_1d_trafo( self,  in, out, degree)
    class(sll_t_maxwell_clamped_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_clamped_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_clamped_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 )
       out(1) = 0._f64
       out(self%n_dofs0) = 0._f64
    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
 
 
  subroutine compute_field_energy( self, efield_dofs1, efield_dofs2, bfield_dofs, energy)
    class(sll_t_maxwell_clamped_1d_trafo)  :: self
    sll_real64, intent( in    )          :: efield_dofs1(:) 
    sll_real64, intent( in    )          :: efield_dofs2(:) 
    sll_real64, intent( in    )          :: bfield_dofs(:) 
    sll_real64, intent(   out )          :: energy 
    !local variables
    sll_real64 :: field_energy(3)
 
 
    field_energy(1) = self%l2norm_squared &
         ( efield_dofs1(1:self%n_dofs1), self%s_deg_1 )
    field_energy(2) = self%l2norm_squared &
         ( efield_dofs2(1:self%n_dofs0), self%s_deg_0 )
    field_energy(3) =self%l2norm_squared &
         ( bfield_dofs(1:self%n_dofs1), self%s_deg_1 )
 
    energy = 0.5_f64*sum(field_energy) 
 
  end subroutine compute_field_energy
 
 
end module sll_m_maxwell_clamped_1d_trafo