sll_m_dg_fields.F90

Go to the documentation of this file.

#define sll_transformation class(sll_c_coordinate_transformation_2d_base)
 
module sll_m_dg_fields
 
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#include "sll_assert.h"
#include "sll_memory.h"
#include "sll_working_precision.h"
 
   use sll_m_ascii_io, only: &
      sll_s_ascii_file_create
 
   use sll_m_cartesian_meshes, only: &
      sll_t_cartesian_mesh_2d
 
   use sll_m_coordinate_transformation_2d_base, only: &
      sll_c_coordinate_transformation_2d_base, &
      sll_p_io_gmsh, &
      sll_p_io_mtv, &
      sll_p_io_xdmf
 
   use sll_m_gauss_lobatto_integration, only: &
      sll_f_gauss_lobatto_points, &
      sll_f_gauss_lobatto_weights
 
   use sll_m_utilities, only: &
      sll_s_int2string
 
   implicit none
 
   public :: &
      sll_s_dg_field_2d_init, &
      sll_t_dg_field_2d, &
      sll_f_new_dg_field_2d
 
   private
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
   type :: sll_t_dg_field_2d
 
      sll_int32                                   :: degree  
      sll_transformation, pointer     :: tau     
      sll_real64, dimension(:, :, :, :), allocatable :: array   
      sll_real64, dimension(:), allocatable :: xgalo   
      sll_real64, dimension(:), allocatable :: wgalo   
      sll_int32                                   :: tag     
      sll_int32                                   :: file_id 
 
   contains
 
      procedure, pass :: init => sll_s_dg_field_2d_init
      procedure, pass :: write_to_file => write_dg_field_2d_to_file
      procedure, pass :: set_value => set_value_dg_field_2d
 
   end type sll_t_dg_field_2d
 
   interface operator(+)
      module procedure dg_field_add
   end interface operator(+)
 
   interface operator(-)
      module procedure dg_field_sub
   end interface operator(-)
 
   sll_int32 :: error
 
contains
 
   subroutine sll_s_dg_field_2d_init(this, degree, tau, init_function)
 
      class(sll_t_dg_field_2d)              :: this
      sll_transformation, target            :: tau           
      sll_real64, external, optional        :: init_function
      sll_int32, intent(in)                 :: degree        
      sll_int32                             :: nc_eta1
      sll_int32                             :: nc_eta2
      sll_int32                             :: error
 
      this%tau => tau
      this%degree = degree
 
      sll_allocate(this%xgalo(degree + 1), error)
      sll_allocate(this%wgalo(degree + 1), error)
 
      this%xgalo = sll_f_gauss_lobatto_points(degree + 1, -1._f64, 1._f64)
      this%wgalo = sll_f_gauss_lobatto_weights(degree + 1)
 
      nc_eta1 = tau%mesh%num_cells1
      nc_eta2 = tau%mesh%num_cells2
      allocate (this%array(1:degree + 1, 1:degree + 1, 1:nc_eta1, 1:nc_eta2))
 
      this%array = 0.0_f64
      if (present(init_function)) then
         call set_value_dg_field_2d(this, init_function, 0.0_f64)
      end if
 
      this%tag = 0
      this%file_id = 0
 
   end subroutine sll_s_dg_field_2d_init
 
   function sll_f_new_dg_field_2d(degree, tau, init_function) result(this)
 
      sll_transformation, pointer           :: tau           
      sll_real64, external, optional        :: init_function
      sll_int32, intent(in)                 :: degree        
      sll_int32                             :: nc_eta1
      sll_int32                             :: nc_eta2
      type(sll_t_dg_field_2d), pointer      :: this
      sll_int32                             :: error
      sll_allocate(this, error)
      this%tau => tau
      this%degree = degree
 
      sll_allocate(this%xgalo(degree + 1), error)
      sll_allocate(this%wgalo(degree + 1), error)
 
      this%xgalo = sll_f_gauss_lobatto_points(degree + 1, -1._f64, 1._f64)
      this%wgalo = sll_f_gauss_lobatto_weights(degree + 1)
 
      associate(lm => this%tau%mesh)
         nc_eta1 = lm%num_cells1
         nc_eta2 = lm%num_cells2
      end associate
 
      allocate (this%array(1:degree + 1, 1:degree + 1, 1:nc_eta1, 1:nc_eta2))
 
      this%array = 0.0_f64
      if (present(init_function)) then
         call set_value_dg_field_2d(this, init_function, 0.0_f64)
      end if
 
      this%tag = 0
      this%file_id = 0
 
   end function sll_f_new_dg_field_2d
 
   subroutine set_value_dg_field_2d(this, init_function, time)
 
      class(sll_t_dg_field_2d) :: this
      sll_real64, external     :: init_function
      sll_real64               :: time
      sll_real64               :: offset(2)
      sll_real64               :: eta1
      sll_real64               :: eta2
      sll_int32                :: i, j, ii, jj
 
      sll_assert(allocated(this%array))
 
      do j = 1, this%tau%mesh%num_cells2
      do i = 1, this%tau%mesh%num_cells1
         offset(1) = this%tau%mesh%eta1_min + real(i - 1, f64)*this%tau%mesh%delta_eta1
         offset(2) = this%tau%mesh%eta2_min + real(j - 1, f64)*this%tau%mesh%delta_eta2
         do jj = 1, this%degree + 1
         do ii = 1, this%degree + 1
            eta1 = offset(1) + 0.5_f64*(this%xgalo(ii) + 1.0_f64)*this%tau%mesh%delta_eta1
            eta2 = offset(2) + 0.5_f64*(this%xgalo(jj) + 1.0_f64)*this%tau%mesh%delta_eta2
            this%array(ii, jj, i, j) = init_function(this%tau%x1(eta1, eta2), &
                                                     this%tau%x2(eta1, eta2), &
                                                     time)
         end do
         end do
      end do
      end do
 
   end subroutine set_value_dg_field_2d
 
   subroutine write_dg_field_2d_to_file(this, field_name, file_format, time)
 
      class(sll_t_dg_field_2d) :: this
      character(len=*)         :: field_name
      sll_int32, optional      :: file_format
      sll_real64, optional     :: time
 
      if (present(file_format)) then
         select case (file_format)
         case (sll_p_io_gmsh)
            call plot_dg_field_2d_with_gmsh(this, field_name)
         case (sll_p_io_mtv)
            call plot_dg_field_2d_with_plotmtv(this, field_name)
         case (sll_p_io_xdmf)
            if (present(time)) then
               call plot_dg_field_2d_with_xdmf(this, field_name, time)
            else
               call plot_dg_field_2d_with_xdmf(this, field_name)
            end if
         case default
            print "(a)", field_name//' write_to_file : Unknown format'
         end select
      else
         call plot_dg_field_2d_with_gnuplot(this, field_name)
      end if
 
   end subroutine write_dg_field_2d_to_file
 
   subroutine plot_dg_field_2d_with_gnuplot(this, field_name)
 
      class(sll_t_dg_field_2d) :: this
      character(len=*)         :: field_name
      sll_int32                :: file_id
      sll_int32                :: gnu_id
      sll_real64               :: eta1, eta2
      sll_real64               :: offset(2)
      sll_int32                :: i, j, ii, jj
      character(len=4)         :: ctag
 
      call sll_s_int2string(this%tag, ctag)
 
      gnu_id = this%file_id
 
      if (gnu_id == 0) then
         call sll_s_ascii_file_create(field_name//".gnu", gnu_id, error)
         rewind(gnu_id)
      else
         open (unit=gnu_id, file=field_name//".gnu", position="append")
      end if
      write (gnu_id, "(a)") "unset key"
      write (gnu_id, "(a)") "set title '"//field_name//" step "//ctag//"'"
      write (gnu_id, "(a)") "splot '"//field_name//ctag//".dat' w l"
 
      call sll_s_ascii_file_create(field_name//ctag//".dat", file_id, error)
 
      associate(lm => this%tau%mesh)
 
         do j = 1, lm%num_cells2
         do i = 1, lm%num_cells1
 
            offset(1) = lm%eta1_min + real(i - 1, f64)*lm%delta_eta1
            offset(2) = lm%eta2_min + real(j - 1, f64)*lm%delta_eta2
            do jj = 1, this%degree + 1
            do ii = 1, this%degree + 1
               eta1 = offset(1) + 0.5_f64*(this%xgalo(ii) + 1.0_f64)*lm%delta_eta1
               eta2 = offset(2) + 0.5_f64*(this%xgalo(jj) + 1.0_f64)*lm%delta_eta2
               write (file_id, *) this%tau%x1(eta1, eta2), &
                  this%tau%x2(eta1, eta2), &
                  sngl(this%array(ii, jj, i, j))
            end do
            write (file_id, *)
            end do
            write (file_id, *)
 
         end do
         end do
 
      end associate
 
      close (file_id)
 
      write (gnu_id, *)
      close (gnu_id)
 
      this%tag = this%tag + 1
      this%file_id = gnu_id
 
   end subroutine plot_dg_field_2d_with_gnuplot
 
   function dg_field_add(W1, W2) result(W3)
 
      type(sll_t_dg_field_2d), intent(in) :: w1
      type(sll_t_dg_field_2d), intent(in) :: w2
      type(sll_t_dg_field_2d)             :: w3
 
      sll_assert(w1%degree == w2%degree)
      sll_assert(allocated(w1%array))
      sll_assert(allocated(w2%array))
 
      w3%array = w1%array + w2%array
 
   end function dg_field_add
 
   function dg_field_sub(W1, W2) result(W3)
 
      type(sll_t_dg_field_2d), intent(in) :: w1
      type(sll_t_dg_field_2d), intent(in) :: w2
      type(sll_t_dg_field_2d)             :: w3
 
      w3%array = w1%array - w2%array
 
   end function dg_field_sub
 
   subroutine plot_dg_field_2d_with_gmsh(this, field_name)
 
      class(sll_t_dg_field_2d) :: this
      character(len=*)         :: field_name
 
      sll_int32                :: file_id
      sll_int32, parameter     :: nlocmax = 16
      sll_int32                :: typelem
      sll_int32, dimension(4)  :: invpermut1 = (/1, 2, 4, 3/)
      sll_int32, dimension(9)  :: invpermut2 = (/1, 5, 2, 8, 9, 6, 4, 7, 3/)
      sll_int32, dimension(16) :: invpermut3 = (/1, 5, 6, 2, 12, 13, 14, 7, 11, 16, 15, 8, 4, 10, 9, 3/)
      sll_int32                :: invpermut(nlocmax), permut(nlocmax)
      sll_int32                :: nloc, nel, neq, iel, ino, inoloc
      sll_int32                :: i, j, k, l, ii, jj, ni, nj
      sll_int32, allocatable  :: connec(:, :)
      sll_real64, allocatable  :: coords(:, :)
      sll_real64, allocatable  :: values(:)
      sll_real64               :: offset(2), eta1, eta2
      character(len=32)        :: my_fmt
 
      if (this%degree > 3) then
         write (*, *) é'ordre non prvu'
         stop
      end if
 
      associate(lm => this%tau%mesh)
 
         ni = lm%num_cells1
         nj = lm%num_cells2
         nloc = (this%degree + 1)*(this%degree + 1)
         nel = ni*nj
         neq = (ni*this%degree + 1)*(nj*this%degree + 1)
 
         sll_allocate(connec(1:nloc, 1:neq), error)
         sll_allocate(coords(1:2, 1:neq), error)
         sll_allocate(values(1:neq), error)
 
         do j = 0, nj - 1
         do i = 0, ni - 1
            iel = 1 + i + j*ni
            offset(1) = lm%eta1_min + i*lm%delta_eta1
            offset(2) = lm%eta2_min + j*lm%delta_eta2
            do jj = 0, this%degree
            do ii = 0, this%degree
               inoloc = 1 + ii + (this%degree + 1)*jj
               ino = 1 + ii + this%degree*i + (this%degree*ni + 1)*(jj + this%degree*j)
               connec(inoloc, iel) = ino
               eta1 = offset(1) + .5_f64*(this%xgalo(ii + 1) + 1._f64)*lm%delta_eta1
               eta2 = offset(2) + .5_f64*(this%xgalo(jj + 1) + 1._f64)*lm%delta_eta2
               coords(1, ino) = this%tau%x1(eta1, eta2)
               coords(2, ino) = this%tau%x2(eta1, eta2)
               values(ino) = this%array(ii + 1, jj + 1, i + 1, j + 1)
            end do
            end do
         end do
         end do
 
         select case (this%degree)
         case (1)
            invpermut(1:nloc) = invpermut1
            typelem = 3
         case (2)
            invpermut(1:nloc) = invpermut2
            typelem = 10
         case (3)
            invpermut(1:nloc) = invpermut3
            typelem = 36
         case default
            write (*, *) é'ordre non prvu'
            stop
         end select
 
         do i = 1, nloc
            permut(invpermut(i)) = i
         end do
 
         call sll_s_ascii_file_create(field_name//".msh", file_id, error)
         write (file_id, '(A)') '$MeshFormat'
         write (file_id, '(A)') '2 0 8'
         write (file_id, '(A)') '$EndMeshFormat'
         write (file_id, '(A)') '$Nodes'
         write (file_id, *) neq
         do k = 1, neq
            write (file_id, '(i6,2x,3f10.5)') k, coords(1:2, k), 0.0
         end do
         write (file_id, '(A)') '$EndNodes'
         write (file_id, '(A)') '$Elements'
         write (file_id, *) nel
 
         write (my_fmt, '(a, i0, a)') '(3i6,', nloc, 'i5)'
 
         k = 0
         do i = 1, lm%num_cells1
         do j = 1, lm%num_cells2
            k = k + 1
            write (file_id, trim(my_fmt)) k, typelem, 0, (connec(l, k), l=1, nloc)
         end do
         end do
         write (file_id, '(A)') '$EndElements'
         write (file_id, '(A)') '$NodeData'
         write (file_id, *) 1
         write (file_id, *) 'values'
         write (file_id, *) 1
         write (file_id, *) 0
         write (file_id, *) 3
         write (file_id, *) 0
         write (file_id, *) 1
         write (file_id, *) neq
         do k = 1, neq
            write (file_id, '(i6,2x,f10.5)') k, values(k)
         end do
         write (file_id, '(A)') '$EndNodeData'
 
         close (file_id)
      end associate
 
   end subroutine plot_dg_field_2d_with_gmsh
 
   subroutine plot_dg_field_2d_with_plotmtv(this, field_name)
 
      class(sll_t_dg_field_2d) :: this
      character(len=*)         :: field_name
 
      sll_int32                :: file_id
      sll_int32                :: ni, nj, ino
      sll_int32                :: i, j, ii, jj
      sll_real64               :: offset(2)
      sll_real64               :: eta1, eta2
 
      offset = 0.0_f64
      call sll_s_ascii_file_create(field_name//".mtv", file_id, error)
 
      write (file_id, *) "$DATA=CONTCURVE"
      write (file_id, *) "%equalscale=T"
      write (file_id, *) "%contfill"
      write (file_id, *) "%toplabel='"//field_name//"'"
 
      associate(lm => this%tau%mesh)
 
         do j = 1, lm%num_cells2
            offset(2) = lm%eta2_min + real(j - 1, f64)*lm%delta_eta2
            do jj = 1, merge(this%degree, this%degree + 1, j < lm%num_cells2)
               do i = 1, lm%num_cells1
                  offset(1) = lm%eta1_min + real(i - 1, f64)*lm%delta_eta1
                  do ii = 1, merge(this%degree, this%degree + 1, i < lm%num_cells1)
                     eta1 = offset(1) + .5_f64*(this%xgalo(ii) + 1._f64)*lm%delta_eta1
                     eta2 = offset(2) + .5_f64*(this%xgalo(jj) + 1._f64)*lm%delta_eta2
                     write (file_id, *) sngl(this%tau%x1(eta1, eta2)), &
                        sngl(this%tau%x2(eta1, eta2)), &
                        sngl(this%array(ii, jj, i, j))
                  end do
               end do
            end do
         end do
 
         write (file_id, *)
 
         write (file_id, *) "$DATA=CURVE3D"
         write (file_id, *) "%equalscale=T"
         write (file_id, *) "%meshplot"
         write (file_id, *) "%toplabel='"//field_name//"'"
 
         do j = 1, lm%num_cells2
         do i = 1, lm%num_cells1
            offset(1) = lm%eta1_min + (i - 1)*lm%delta_eta1
            offset(2) = lm%eta2_min + (j - 1)*lm%delta_eta2
            do jj = 1, this%degree
            do ii = 1, this%degree
               eta1 = offset(1) + .5_f64*(this%xgalo(ii) + 1._f64)*lm%delta_eta1
               eta2 = offset(2) + .5_f64*(this%xgalo(jj) + 1._f64)*lm%delta_eta2
               write (file_id, *) sngl(this%tau%x1(eta1, eta2)), &
                  sngl(this%tau%x2(eta1, eta2)), &
                  sngl(this%array(ii, jj, i, j))
               eta1 = offset(1) + .5_f64*(this%xgalo(ii + 1) + 1._f64)*lm%delta_eta1
               eta2 = offset(2) + .5_f64*(this%xgalo(jj) + 1._f64)*lm%delta_eta2
               write (file_id, *) sngl(this%tau%x1(eta1, eta2)), &
                  sngl(this%tau%x2(eta1, eta2)), &
                  sngl(this%array(ii + 1, jj, i, j))
               eta1 = offset(1) + .5_f64*(this%xgalo(ii + 1) + 1._f64)*lm%delta_eta1
               eta2 = offset(2) + .5_f64*(this%xgalo(jj + 1) + 1._f64)*lm%delta_eta2
               write (file_id, *) sngl(this%tau%x1(eta1, eta2)), &
                  sngl(this%tau%x2(eta1, eta2)), &
                  sngl(this%array(ii + 1, jj + 1, i, j))
               eta1 = offset(1) + .5_f64*(this%xgalo(ii) + 1._f64)*lm%delta_eta1
               eta2 = offset(2) + .5_f64*(this%xgalo(jj + 1) + 1._f64)*lm%delta_eta2
               write (file_id, *) sngl(this%tau%x1(eta1, eta2)), &
                  sngl(this%tau%x2(eta1, eta2)), &
                  sngl(this%array(ii, jj + 1, i, j))
               eta1 = offset(1) + .5_f64*(this%xgalo(ii) + 1._f64)*lm%delta_eta1
               eta2 = offset(2) + .5_f64*(this%xgalo(jj) + 1._f64)*lm%delta_eta2
               write (file_id, *) sngl(this%tau%x1(eta1, eta2)), &
                  sngl(this%tau%x2(eta1, eta2)), &
                  sngl(this%array(ii, jj, i, j))
               write (file_id, *)
            end do
            end do
         end do
         end do
 
         ni = lm%num_cells1
         nj = lm%num_cells2
 
         do j = 0, nj - 1
         do i = 0, ni - 1
            offset(1) = lm%eta1_min + i*lm%delta_eta1
            offset(2) = lm%eta2_min + j*lm%delta_eta2
            do jj = 0, this%degree
            do ii = 0, this%degree
               ino = 1 + ii + this%degree*i + (this%degree*ni + 1)*(jj + this%degree*j)
               eta1 = offset(1) + .5*(this%xgalo(ii + 1) + 1.)*lm%delta_eta1
               eta2 = offset(2) + .5*(this%xgalo(jj + 1) + 1.)*lm%delta_eta2
               write (file_id, "(a)", advance="no") "@text x1="
               write (file_id, "(g15.3)", advance="no") this%tau%x1(eta1, eta2)
               write (file_id, "(a)", advance="no") " y1="
               write (file_id, "(g15.3)", advance="no") this%tau%x2(eta1, eta2)
               write (file_id, "(a)", advance="no") " z1=0. lc=5 ll='"
               write (file_id, "(i4)", advance="no") ino
               write (file_id, "(a)") "'"
            end do
            end do
         end do
         end do
 
         write (file_id, *) "$end"
         close (file_id)
      end associate
 
   end subroutine plot_dg_field_2d_with_plotmtv
 
   subroutine plot_dg_field_2d_with_xdmf(this, field_name, time)
 
      class(sll_t_dg_field_2d) :: this
      character(len=*)         :: field_name
      sll_real64, optional     :: time
 
      sll_int32                :: file_id
      sll_int32                :: i, j, k, ii, jj
      sll_real64               :: offset(2)
      sll_real64               :: eta1, eta2
      sll_int32                :: clength
 
      clength = len_trim(field_name)
 
      sll_assert(this%degree < 10)
 
      call sll_s_ascii_file_create(field_name//".xmf", file_id, error)
 
      write (file_id, "(a)") "<?xml version='1.0' ?>"
      write (file_id, "(a)") "<!DOCTYPE Xdmf SYSTEM 'Xdmf.dtd' []>"
      write (file_id, "(a)") "<Xdmf Version='2.'>"
      write (file_id, "(a)") "<Domain>"
      write (file_id, "(a)") "<Grid Name='Mesh' GridType='Collection'>"
      if (present(time)) then
         write (file_id, "(a13,g15.3,a3)") "<Time Value='", time, "'/>"
      end if
 
      associate(lm => this%tau%mesh)
 
         k = 0
         do i = 1, lm%num_cells1
         do j = 1, lm%num_cells2
            offset(1) = lm%eta1_min + real(i - 1, f64)*lm%delta_eta1
            offset(2) = lm%eta2_min + real(j - 1, f64)*lm%delta_eta2
            k = k + 1
            write (file_id, "(a,i6,a)") "<Grid Name='Mesh", k, "' GridType='Uniform'>"
            write (file_id, "(a,2i5,a)") &
               "<Topology TopologyType='2DSMesh' NumberOfElements='", &
               this%degree + 1, this%degree + 1, "'/>"
            write (file_id, "(a)") "<Geometry GeometryType='X_Y'>"
            write (file_id, "(a,2i3,a)") "<DataItem Dimensions='", this%degree + 1, &
               this%degree + 1, "' NumberType='Float' Format='XML'>"
            do jj = 1, this%degree + 1
            do ii = 1, this%degree + 1
               eta1 = offset(1) + .5*(this%xgalo(ii) + 1.0_f64)*lm%delta_eta1
               eta2 = offset(2) + .5*(this%xgalo(jj) + 1.0_f64)*lm%delta_eta2
               write (file_id, "(f7.3)", advance='no') sngl(this%tau%x1(eta1, eta2))
            end do
            write (file_id, *)
            end do
            write (file_id, "(a)") "</DataItem>"
            write (file_id, "(a,2i3,a)") "<DataItem Dimensions='", this%degree + 1, &
               this%degree + 1, "' NumberType='Float' Format='XML'>"
            do jj = 1, this%degree + 1
            do ii = 1, this%degree + 1
               eta1 = offset(1) + .5*(this%xgalo(ii) + 1.0_f64)*lm%delta_eta1
               eta2 = offset(2) + .5*(this%xgalo(jj) + 1.0_f64)*lm%delta_eta2
               write (file_id, "(f7.3)", advance='no') sngl(this%tau%x2(eta1, eta2))
            end do
            write (file_id, *)
            end do
            write (file_id, "(a)") "</DataItem>"
            write (file_id, "(a)") "</Geometry>"
            write (file_id, "(a)") "<Attribute Name='values' &
                                 & AttributeType='Scalar' Center='Node'>"
            write (file_id, "(a,2i3,a)") "<DataItem Dimensions='", &
               this%degree + 1, this%degree + 1, &
               "' NumberType='Float' Precision='4' Format='XML'>"
            do jj = 1, this%degree + 1
            do ii = 1, this%degree + 1
               eta1 = offset(1) + .5_f64*(this%xgalo(ii) + 1._f64)*lm%delta_eta1
               eta2 = offset(2) + .5_f64*(this%xgalo(jj) + 1._f64)*lm%delta_eta2
               write (file_id, "(f7.3)", advance='no') sngl(this%array(ii, jj, i, j))
            end do
            write (file_id, *)
            end do
            write (file_id, "(a)") "</DataItem>"
            write (file_id, "(a)") "</Attribute>"
            write (file_id, "(a)") "</Grid>"
         end do
         end do
         write (file_id, "(a)") "</Grid>"
         write (file_id, "(a)") "</Domain>"
         write (file_id, "(a)") "</Xdmf>"
 
      end associate
 
   end subroutine plot_dg_field_2d_with_xdmf
 
end module sll_m_dg_fields