sll_m_advection_1d_PSM.F90

Go to the documentation of this file.

!**************************************************************
!  Copyright INRIA
!  Authors :
!     CALVI project team
!
!  This code SeLaLib (for Semi-Lagrangian-Library)
!  is a parallel library for simulating the plasma turbulence
!  in a tokamak.
!
!  This software is governed by the CeCILL-B license
!  under French law and abiding by the rules of distribution
!  of free software.  You can  use, modify and redistribute
!  the software under the terms of the CeCILL-B license as
!  circulated by CEA, CNRS and INRIA at the following URL
!  "http://www.cecill.info".
!**************************************************************
 
! in development; should be at least cubic splines
! attached with computation of characteristics
 
module sll_m_advection_1d_psm
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#include "sll_assert.h"
#include "sll_memory.h"
#include "sll_working_precision.h"
 
   use sll_m_advection_1d_base, only: &
      sll_c_advector_1d
 
   implicit none
 
   public :: &
      sll_f_new_psm_1d_advector
 
   private
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
   type, extends(sll_c_advector_1d) :: psm_1d_advector
      sll_real64, dimension(:), pointer :: buf1d
      sll_real64, dimension(:), pointer :: buf1d_out
      sll_real64, dimension(:), pointer :: dtab
      sll_real64, dimension(:), pointer :: ltab
      sll_real64, dimension(:), pointer :: mtab
      sll_real64, dimension(:), pointer :: alphax
      sll_real64, dimension(:), pointer :: prim
      sll_real64, dimension(:), pointer :: bufout0
      sll_real64, dimension(:), pointer :: bufout
      sll_real64, dimension(:), pointer :: p
      sll_real64 :: eta_min
      sll_real64 :: eta_max
      sll_int32 :: npts
      sll_int32 :: nbdr
   contains
      procedure, pass(adv) :: initialize => &
         initialize_psm_1d_advector
      procedure, pass(adv) :: advect_1d => &
         psm_advect_1d
      procedure, pass(adv) :: advect_1d_constant => &
         psm_advect_1d_constant
      procedure, pass(adv) :: delete => delete_psm_1d_advector
   end type psm_1d_advector
 
contains
   function sll_f_new_psm_1d_advector( &
      Npts, &
      eta_min, &
      eta_max, &
      Nbdr) &
      result(adv)
      type(psm_1d_advector), pointer :: adv
      sll_int32, intent(in) :: npts
      sll_real64, intent(in) :: eta_min
      sll_real64, intent(in) :: eta_max
      sll_int32, intent(in), optional :: nbdr
      sll_int32 :: ierr
 
      sll_allocate(adv, ierr)
 
      call initialize_psm_1d_advector( &
         adv, &
         npts, &
         eta_min, &
         eta_max, &
         nbdr)
 
   end function sll_f_new_psm_1d_advector
 
   subroutine initialize_psm_1d_advector( &
      adv, &
      Npts, &
      eta_min, &
      eta_max, &
      Nbdr)
      class(psm_1d_advector), intent(inout) :: adv
      sll_int32, intent(in) :: npts
      sll_real64, intent(in) :: eta_min
      sll_real64, intent(in) :: eta_max
      sll_int32, intent(in), optional :: nbdr
      sll_int32 :: ierr
      sll_int32 :: i
      sll_int32 :: nx
 
      adv%Npts = npts
      adv%eta_min = eta_min
      adv%eta_max = eta_max
 
      sll_allocate(adv%buf1d(npts), ierr)
      sll_allocate(adv%buf1d_out(npts), ierr)
      nx = npts - 1
 
      if (present(nbdr)) then
         adv%Nbdr = nbdr
      else
         adv%Nbdr = 100
      end if
 
      !print *,Nx,adv%Nbdr
 
      sll_allocate(adv%dtab(0:nx - 1), ierr)
      sll_allocate(adv%ltab(0:nx - 2), ierr)
      sll_allocate(adv%mtab(0:nx - 1), ierr)
      sll_allocate(adv%alphax(0:nx - 1), ierr)
      sll_allocate(adv%prim(-adv%Nbdr:nx - 1 + adv%Nbdr), ierr)
      sll_allocate(adv%bufout0(0:nx - 1 + 2*adv%Nbdr), ierr)
      sll_allocate(adv%bufout(0:nx), ierr)
      sll_allocate(adv%p(0:nx - 1), ierr)
 
      adv%dtab(0) = 4._f64
      adv%mtab(0) = 1._f64
      do i = 0, nx - 3
         adv%ltab(i) = 1._f64/adv%dtab(i)
         adv%dtab(i + 1) = 4._f64 - adv%ltab(i)
         adv%mtab(i + 1) = -adv%mtab(i)/adv%dtab(i)
      end do
      adv%ltab(nx - 2) = (1._f64 + adv%mtab(nx - 2))/adv%dtab(nx - 2)
      adv%dtab(nx - 1) = 4._f64 - adv%ltab(nx - 2)*(adv%mtab(nx - 2) + 1._f64)
      do i = 0, nx - 3
         adv%dtab(nx - 1) = adv%dtab(nx - 1) + adv%mtab(i)*adv%mtab(i + 1)
      end do
      do i = 0, nx - 1
         adv%dtab(i) = 1._f64/adv%dtab(i)
      end do
 
!  dtab= (double*) malloc(sizeof(double)*(Nx));ltab= (double*) malloc(sizeof(double)*(Nx-1));
!mtab= (double*) malloc(sizeof(double)*(Nx));
!  dtab[0]=4.;mtab[0]=1.;for(i=0;i<Nx-2 ;i++){ltab[i]=1./dtab[i];dtab[i+1]=4.-ltab[i];mtab[i+1]=-mtab[i]/dtab[i];}
!  ltab[Nx-2]=(1.+mtab[Nx-2])/dtab[Nx-2];dtab[Nx-1]=4-ltab[Nx-2]*(mtab[Nx-2]+1);
!  for(i=0;i<Nx-2;i++)dtab[Nx-1]=dtab[Nx-1]+mtab[i]*mtab[i+1];for(i=0;i<Nx;i++)dtab[i]=1./dtab[i];
 
   end subroutine initialize_psm_1d_advector
 
   subroutine psm_advect_1d( &
      adv, &
      A, &
      dt, &
      input, &
      output)
      class(psm_1d_advector) :: adv
      sll_real64, dimension(:), intent(in) :: a
      sll_real64, intent(in) :: dt
      sll_real64, dimension(:), intent(in) :: input
      sll_real64, dimension(:), intent(out) :: output
      sll_int32 :: npts
      sll_int32 :: nbdr
      sll_int32 :: i
      sll_int32 :: nx
      sll_int32 :: s
      sll_real64 :: dtmp2
      sll_real64 :: x
      sll_real64 :: dx
      sll_int32 :: ix
      sll_int32 :: ix1
      sll_real64 :: result
      sll_real64 :: m
      sll_real64 :: w(0:2)
      sll_real64 :: tmp
      sll_real64 :: df0
      sll_real64 :: df1
      sll_real64 :: dt_loc
 
      dt_loc = 0.5_f64*dt
 
      npts = adv%Npts
 
      nbdr = adv%Nbdr
 
      nx = npts - 1
 
      dtmp2 = 1._f64/(adv%eta_max - adv%eta_min)
 
      dx = 1._f64/real(nx, f64)
 
      adv%alphax(0:nx - 1) = 0._f64
 
      !print *,size(A),Nx,Npts,size(adv%alphax)
 
      do i = 0, nx - 1
         do s = 0, 9
            x = real(i, f64)*dx - adv%alphax(i)
            if (x >= 1._f64) then
               x = x - 1._f64
            end if
            if (x < 0._f64) then
               x = x + 1._f64
            end if
            if (x == 1) then
               x = 0._f64
            end if
            ix = floor(x*real(nx, f64))
            if ((x >= 1._f64) .or. (x < 0._f64)) then
               print *, '#x too big or too small'
               print *, '#in  PSM_advect_1d'
               print *, x, nx, ix, i, s, adv%alphax(i)
               stop
            end if
            ix1 = ix + 1
            if (ix1 == nx) then
               ix1 = 0
            end if
 
            !print *,'ix=',i,s,ix,ix1,x,dtmp2,adv%alphax(i),dt
 
            result = (1._f64 - x)*a(ix + 1) + x*a(ix1 + 1)
            adv%alphax(i) = result*dtmp2*dt_loc
         end do
      end do
 
      dtmp2 = adv%alphax(nx - 1)
      do i = nx - 1, 1, -1
         adv%alphax(i) = adv%alphax(i) + adv%alphax(i - 1)
      end do
      adv%alphax(0) = adv%alphax(0) + dtmp2
 
!  dtmp2=1./(xmax-xmin);dtmp3=1./(ymax-ymin);
!  #ifdef PTFX1D
!  for(i=0;i<Nx*Ny;i++)alphax[i]=0.;
!  //for(s=0;s<10;s++)for(i=0;i<Nx;i++)for(j=1;j<Ny;j++){
!  for(j=0;j<Ny;j++){p=&alphax[Nx*j];for(i=0;i<Nx;i++)for(s=0;s<10;s++){
!    //localization of x,y ix and iy
!    x=i*dx-p[i];if(x>=1.)x-=1.;if(x<0.)x+=1.;if(x>=1.)x-=1.;ix=(int)(x*Nx);
!    if(x>=1. || x<0.){fprintf(stderr,"x too big/small %1.1000lg i=%d j=%d\n",x,i,j);exit(1);}
!    assert(x>=0.&&x<1.);assert(ix>=0 &&ix<Nx);x=x*Nx-ix;assert(x>=0 &&x<1.);
!    ix1=ix+1;if(ix1==Nx)ix1=0;
!    //result=(1.-x)*Ey[j+(Ny+1)*ix]+x*Ey[j+(Ny+1)*ix1];alphax[(i+Nx*j)]=result*dtmp2*dt;
!    result=(1.-x)*Ey[ix+Nx*j]+x*Ey[ix1+Nx*j];p[i]=result*dtmp2*dt;
!  }}
 
!  if(interp==5 || interp==11 || interp==6 || interp<0){
!    for(j=0;j<Ny;j++){p=&alphax[Nx*j];dtmp2=p[Nx-1];for(i=Nx-1;i>0;i--)p[i]=(p[i]+p[i-1]);p[0]=(p[0]+dtmp2);
!  }
 
      adv%p(0:nx - 1) = input(1:nx)
      m = sum(adv%p(0:nx - 1))
      adv%prim(0) = adv%p(0)
      do i = 1, nx - 1
         adv%prim(i) = adv%prim(i - 1) + adv%p(i)
      end do
 
      do i = -nbdr, -1
         adv%prim(i) = adv%prim(i + nx) - m
      end do
      do i = nx, nx + nbdr - 1
         adv%prim(i) = adv%prim(i - nx) + m
      end do
 
      do i = 0, nx - 1
         adv%bufout0(i + nbdr) = 3._f64*(adv%p(i) + adv%p(mod(i + nx - 1, nx)))
      end do
      do i = 1, nx - 1
         adv%bufout0(i + nbdr) = adv%bufout0(i + nbdr) - adv%bufout0(i - 1 + nbdr)*adv%ltab(i - 1)
      end do
      do i = 0, nx - 3
         adv%bufout0(nx - 1 + nbdr) = adv%bufout0(nx - 1 + nbdr) + adv%mtab(i + 1)*adv%bufout0(i + nbdr)
      end do
      adv%bufout0(nx - 1 + nbdr) = adv%bufout0(nx - 1 + nbdr)*adv%dtab(nx - 1)
      adv%bufout0(nx - 2 + nbdr) = adv%dtab(nx - 2)*(adv%bufout0(nx - 2 + nbdr) &
                                                     - (1._f64 + adv%mtab(nx - 2))*adv%bufout0(nx - 1 + nbdr))
      do i = nx - 3, 0, -1
         adv%bufout0(i + nbdr) = adv%dtab(i) &
                                 *(adv%bufout0(i + nbdr) - adv%bufout0(i + 1 + nbdr) - adv%mtab(i)*adv%bufout0(nx - 1 + nbdr))
      end do
      do i = 0, nbdr - 1
         adv%bufout0(i) = adv%bufout0(nx + i)
      end do
      do i = 0, nbdr - 1
         adv%bufout0(nx + nbdr + i) = adv%bufout0(nbdr + i)
      end do
 
!    M=0.;for(i=0;i<Nx;i++)M+=p[i];prim[0]=p[0];
!    for(i=1;i<Nx;i++)prim[i]=prim[i-1]+p[i];
!    for(i=-Nbdr;i<0;i++)prim[i]=prim[i+Nx]-M;for(i=Nx;i<Nx+Nbdr;i++)prim[i]=prim[i-Nx]+M;
!    for(i=0;i<Nx;i++)bufout0[i+Nbdr]=3*(p[i]+p[(i+Nx-1)%Nx]);
!     /*resolution of L phi = ftab*/
!    for(i=1;i<Nx;i++) bufout0[i+Nbdr]-=bufout0[i-1+Nbdr]*ltab[i-1];
!    for(i=0;i<Nx-2;i++)bufout0[Nx-1+Nbdr]+=mtab[i+1]*bufout0[i+Nbdr];
!    /*resolution of U eta =phi*/
!    bufout0[Nx-1+Nbdr]*=dtab[Nx-1];
!    bufout0[Nx-2+Nbdr]=dtab[Nx-2]*(bufout0[Nx-2+Nbdr]-(1.+mtab[Nx-2])*bufout0[Nx-1+Nbdr]);
!    for(i=Nx-3;i>=0;i--) bufout0[i+Nbdr]=dtab[i]*(bufout0[i+Nbdr]-bufout0[i+1+Nbdr]-mtab[i]*bufout0[Nx-1+Nbdr]);
!    for(i=0;i<Nbdr;i++)bufout0[i]=bufout0[Nx+i];
!    for(i=0;i<Nbdr;i++)bufout0[Nx+Nbdr+i]=bufout0[Nbdr+i];
 
      do i = 0, nx
         x = real(i, f64)*dx - adv%alphax(mod(i, nx))
         ix = floor(x*real(nx, f64))
         x = x*real(nx, f64) - real(ix, f64)
         if (x == 1) then
            x = 0._f64
            ix = ix - 1
         end if
         if ((x < 0._f64) .or. (x >= 1._f64)) then
            print *, '#Localization problem in PSM_advect_1d'
            print *, x, ix, nx, adv%alphax(mod(i, nx))
            stop
         end if
         df0 = adv%bufout0(ix + nbdr)
         df1 = adv%bufout0(ix + 1 + nbdr)
         tmp = adv%p(mod(ix + nx, nx))
         w(0) = x*(x - 1._f64)*(x - 1._f64)
         w(1) = x*x*(x - 1._f64)
         w(2) = x*x*(3._f64 - 2._f64*x)
         result = w(0)*df0 + w(1)*df1 + w(2)*tmp
         adv%bufout(i) = result + adv%prim(ix - 1)
      end do
 
      do i = 0, nx - 1
         output(i + 1) = adv%bufout(i + 1) - adv%bufout(i)
      end do
 
      output(nx + 1) = output(1)
 
!    for(i=0;i<Nx+1;i++){
!      //i1=(i-1+Nx)%Nx;
!      x=(i)*dx-alphax[i%Nx];
!      ix=floor(x*Nx);
!      assert(x>-1.&&x<2.);assert(ix-1>=-Nbdr &&ix+1<Nx+Nbdr);x=x*Nx-ix;assert(x>=0 &&x<1.);
!      if(x>=2. || x<=-1.){fprintf(stderr,"x too big/small %1.1000lg i=%d j=%d\n",x,i,j);exit(1);}
!      df0=bufout0[ix+Nbdr];df1=bufout0[ix+1+Nbdr];tmp=p[(ix+Nx)%Nx];
!      lim(tmp,df0,df1);
!      w[0]=x*(x-1)*(x-1);
!      w[1]=x*x*(x-1);
!      w[2]=x*x*(3-2*x);
!      result=w[0]*df0+w[1]*df1+w[2]*tmp;
!      bufout[i]=result+prim[ix-1];
!    }
!for(i=0;i<Nx;i++)p[i]=bufout[i+1]-bufout[i];
 
   end subroutine psm_advect_1d
 
   subroutine psm_advect_1d_constant( &
      adv, &
      A, &
      dt, &
      input, &
      output)
      class(psm_1d_advector) :: adv
      sll_real64, intent(in) :: a
      sll_real64, intent(in) :: dt
      sll_real64, dimension(:), intent(in) :: input
      sll_real64, dimension(:), intent(out) :: output
 
      !this version is not optimized
 
      return
      sll_assert(storage_size(adv) > 0)
      output = input + a*dt
 
   end subroutine psm_advect_1d_constant
 
   subroutine delete_psm_1d_advector(adv)
      class(psm_1d_advector), intent(inout) :: adv
      sll_int32 :: ierr
      sll_deallocate(adv%buf1d, ierr)
      sll_deallocate(adv%buf1d_out, ierr)
      sll_deallocate(adv%dtab, ierr)
      sll_deallocate(adv%ltab, ierr)
      sll_deallocate(adv%mtab, ierr)
      sll_deallocate(adv%alphax, ierr)
      sll_deallocate(adv%prim, ierr)
      sll_deallocate(adv%bufout0, ierr)
      sll_deallocate(adv%bufout, ierr)
      sll_deallocate(adv%p, ierr)
   end subroutine delete_psm_1d_advector
 
end module sll_m_advection_1d_psm