mod_initialize.f90 Source File

!> This module handles the initialization of various components of amrvac
module mod_initialize
  use mod_comm_lib, only: mpistop

  implicit none
  private

  logical :: initialized_already = .false.

  ! Public methods
  public :: initialize_amrvac

contains

  !> Initialize amrvac: read par files and initialize variables
  subroutine initialize_amrvac()
    use mod_global_parameters
    use mod_input_output
    use mod_physics, only: phys_check, phys_check_params
    use mod_usr_methods, only: usr_set_parameters
    use mod_bc_data, only: bc_data_init
    use mod_init_datafromfile, only: read_data_init
    use mod_comm_lib, only: init_comm_types
    use mod_connectivity, only: nbprocs_info

    if (initialized_already) return

    ! Check whether the user has loaded a physics module
    call phys_check()

    ! Read input files
    call read_par_files()
    call initialize_vars()
    call init_comm_types()
    call nbprocs_info%init(npe=npe)

    ! Possibly load boundary condition data or initial data
    call bc_data_init()
    call read_data_init()

    if(associated(usr_set_parameters)) call usr_set_parameters()

    call phys_check_params()

    initialized_already = .true.
  end subroutine initialize_amrvac

  !> Initialize (and allocate) simulation and grid variables
  subroutine initialize_vars
    use mod_forest
    use mod_global_parameters
    use mod_ghostcells_update
    use mod_fix_conserve, only: pflux
    use mod_amr_fct, only: pface, fine_neighbors, old_neighbor
    use mod_geometry

    integer :: igrid, level, ipe, ig1,ig2,ig3
    integer :: nx1, nx2, nx3, iside 
    logical :: ok

    allocate(ps(max_blocks))
    allocate(ps1(max_blocks))
    allocate(ps2(max_blocks))
    allocate(ps3(max_blocks))
    allocate(ps4(max_blocks))
    
    allocate(psc(max_blocks))
    
    
    allocate(ps_sub(max_blocks+1)) ! since we reserve one for communication
    allocate(neighbor(2,-1:1,-1:1,-1:1,max_blocks),neighbor_child(2,0:3,0:3,&
       0:3,max_blocks))
    allocate(neighbor_type(-1:1,-1:1,-1:1,max_blocks),neighbor_active(-1:1,&
       -1:1,-1:1,max_blocks))
    allocate(neighbor_pole(-1:1,-1:1,-1:1,max_blocks))
    allocate(igrids(max_blocks),igrids_active(max_blocks),&
       igrids_passive(max_blocks), idphyb(ndim,max_blocks) )
    allocate(rnode(rnodehi,max_blocks),rnode_sub(rnodehi,max_blocks))
    allocate(node(nodehi,max_blocks),node_sub(nodehi,max_blocks),&
       phyboundblock(max_blocks))

    allocate( bg(1:nstep) )
    
    do istep = 1 , nstep
       bg(istep)%istep = istep
       allocate( bg(istep)%w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3, 1:nw,&
            1:max_blocks) )
       
       
    end do
    allocate( bgc(1) )
    ixCoGmin1=1;ixCoGmin2=1;ixCoGmin3=1;
    ixCoGmax1=(ixGhi1-2*nghostcells)/2+2*nghostcells
    ixCoGmax2=(ixGhi2-2*nghostcells)/2+2*nghostcells
    ixCoGmax3=(ixGhi3-2*nghostcells)/2+2*nghostcells;

    allocate ( bgc(1)%w(ixCoGmin1:ixCoGmax1, ixCoGmin2:ixCoGmax2, &
    ixCoGmin3:ixCoGmax3, 1:nw, 1:max_blocks) )
     
     

    ! The cell metrics follow the same layout as bg%w: one allocation for all
    ! blocks with the grid index last. They are built per block on the device
    ! by fill_geometry_device, which assumes a uniform block - it derives every
    ! metric analytically from the block's corner and spacing in rnode. The one
    ! non-uniform radial mesh supported, LOG_RADIUS, is still uniform in the
    ! logical coordinate xi = ln(r + log_r0), so the assumption holds; general
    ! grid stretching was removed rather than given a device path.
    call set_geometry_index_ranges()
    call alloc_geometry(bgeo, ixGlo1,ixGlo2,ixGlo3, ixGhi1,ixGhi2,ixGhi3, &
         ixGextmin1,ixGextmin2,ixGextmin3, ixGextmax1,ixGextmax2,ixGextmax3)
    call alloc_geometry(bgeoc, ixCoGmin1,ixCoGmin2,ixCoGmin3, ixCoGmax1,&
         ixCoGmax2,ixCoGmax3, ixCoGmin1,ixCoGmin2,ixCoGmin3, ixCoGmax1,&
         ixCoGmax2,ixCoGmax3)
    
    ! What exists is device-resident, because the device is what produces it:
    ! fill_geometry_device builds every allocated member of geo_t there, dx
    ! included even though no kernel reads it - it has to live where it is
    ! written, and is pulled back only on demand.
    
    

    do igrid = 1, max_blocks
       ps(igrid)%igrid  = igrid; ps(igrid)%istep  = 1
       ps1(igrid)%igrid = igrid; ps1(igrid)%istep = 2
       ps2(igrid)%igrid = igrid; ps2(igrid)%istep = 3
       ps3(igrid)%igrid = igrid; ps3(igrid)%istep = 4
       ps4(igrid)%igrid = igrid; ps4(igrid)%istep = 5
    end do

    
    ! allocate mesh for particles
    if(use_particles) allocate(gridvars(max_blocks))
    if(stagger_grid) then
      allocate(pface(2,3,max_blocks),fine_neighbors(2,2,2,3,max_blocks))
      allocate(old_neighbor(2,-1:11,-1:12,-1:13,max_blocks))
    end if

    it=it_init
    global_time=time_init

    dt=zero

    ! no poles initially
    neighbor_pole=0

    ! check resolution
    if (mod(ixGhi1,2)/=0.or.mod(ixGhi2,2)/=0.or.mod(ixGhi3,2)/=0) then
       call mpistop("mesh widths must give even number grid points")
    end if
    ixMlo1=ixGlo1+nghostcells;ixMlo2=ixGlo2+nghostcells
    ixMlo3=ixGlo3+nghostcells;ixMhi1=ixGhi1-nghostcells
    ixMhi2=ixGhi2-nghostcells;ixMhi3=ixGhi3-nghostcells;
    

    nx1 = ixMhi1-ixMlo1+1
    nx2 = ixMhi2-ixMlo2+1
    nx3 = ixMhi3-ixMlo3+1

    allocate(pflux(2,3))
    
    
    do iside = 1, 2
       
      allocate(pflux(iside,1)%flux(1,1:nx2,1:nx3,1:nwflux,1:max_blocks))
      pflux(iside,1)%flux = 0.0d0
      
         
      allocate(pflux(iside,2)%flux(1:nx1,1,1:nx3,1:nwflux,1:max_blocks))
      pflux(iside,2)%flux = 0.0d0
      
         
      allocate(pflux(iside,3)%flux(1:nx1,1:nx2,1,1:nwflux,1:max_blocks))
      pflux(iside,3)%flux = 0.0d0
      

    end do

    if (nbufferx1>(ixMhi1-ixMlo1+1).or.nbufferx2>(ixMhi2-ixMlo2+&
       1).or.nbufferx3>(ixMhi3-ixMlo3+1)) then
       write(unitterm,*) "nbufferx^D bigger than mesh size makes no sense."
       write(unitterm,*) "Decrease nbufferx or increase mesh size"
       call mpistop("")
    end if

    ! initialize dx arrays on finer (>1) levels
    do level=2,refine_max_level
       dx(1,level) = dx(1,level-1) * half
       dx(2,level) = dx(2,level-1) * half
       dx(3,level) = dx(3,level-1) * half  ! refine ratio 2
    end do
    
    ! domain decomposition
    ! physical extent of a grid block at level 1, per dimension
    dg1(1)=dx(1,1)*dble(block_nx1)
    dg2(1)=dx(2,1)*dble(block_nx2)
    dg3(1)=dx(3,1)*dble(block_nx3)
    ! number of grid blocks at level 1 in simulation domain, per dimension
    ng1(1)=nint((xprobmax1-xprobmin1)/dg1(1))
    ng2(1)=nint((xprobmax2-xprobmin2)/dg2(1))
    ng3(1)=nint((xprobmax3-xprobmin3)/dg3(1))
    ! total number of grid blocks at level 1
    nglev1=ng1(1)*ng2(1)*ng3(1)

    do level=2,refine_max_level
       dg1(level)=half*dg1(level-1);dg2(level)=half*dg2(level-1)
       dg3(level)=half*dg3(level-1);
       ng1(level)=ng1(level-1)*2;ng2(level)=ng2(level-1)*2
       ng3(level)=ng3(level-1)*2;
    end do

    ! check that specified stepsize correctly divides domain
    ok=((abs(dble(ng1(1))*dg1(1)-(xprobmax1-xprobmin1))<=&
       smalldouble).and.(abs(dble(ng2(1))*dg2(1)-(xprobmax2-xprobmin2))<=&
       smalldouble).and.(abs(dble(ng3(1))*dg3(1)-(xprobmax3-xprobmin3))<=&
       smalldouble))
    if (.not.ok) then
       write(unitterm,*)"domain cannot be divided by meshes of given gridsize"
       call mpistop("domain cannot be divided by meshes of given gridsize")
    end if

    poleB=.false.
    if (.not.slab) call set_pole
    call check_pole_setup
    ! poleB is declare create'd but had no device copy; nothing on the device
    ! reads it today, and an uninitialised one is a trap for whoever first does
    

    ! number of grid blocks at level 1 along a dimension, which does not have a pole or periodic boundary, 
    ! must be larger than 1 for a rectangular AMR mesh
    if((ng1(1)/=1.or.ng2(1)/=1.or.ng3(1)/=1).and.refine_max_level>1) then
      
      if(ng1(1)==1.and..not.poleB(1,1).and..not.poleB(2,&
         1).and..not.periodB(1).and..not.aperiodB(1)) then
        write(unitterm,"(a,i2,a)")&
            "number of grid blocks at level 1 in dimension",1,&
           " be larger than 1 for a rectangular AMR mesh!"
        write(unitterm,"(a,i1)") "increase domain_nx",1
        call mpistop("")
      end if
      
      
      if(ng2(1)==1.and..not.poleB(1,2).and..not.poleB(2,&
         2).and..not.periodB(2).and..not.aperiodB(2)) then
        write(unitterm,"(a,i2,a)")&
            "number of grid blocks at level 1 in dimension",2,&
           " be larger than 1 for a rectangular AMR mesh!"
        write(unitterm,"(a,i1)") "increase domain_nx",2
        call mpistop("")
      end if
      
      
      if(ng3(1)==1.and..not.poleB(1,3).and..not.poleB(2,&
         3).and..not.periodB(3).and..not.aperiodB(3)) then
        write(unitterm,"(a,i2,a)")&
            "number of grid blocks at level 1 in dimension",3,&
           " be larger than 1 for a rectangular AMR mesh!"
        write(unitterm,"(a,i1)") "increase domain_nx",3
        call mpistop("")
      end if
      
    end if

    ! initialize connectivity data
    igridstail=0

    ! allocate memory for forest data structures
    allocate(level_head(refine_max_level),level_tail(refine_max_level))
    do level=1,refine_max_level
       nullify(level_head(level)%node,level_tail(level)%node)
    end do

    allocate(igrid_to_node(max_blocks,0:npe-1))
    do ipe=0,npe-1
       do igrid=1,max_blocks
          nullify(igrid_to_node(igrid,ipe)%node)
       end do
    end do

    allocate(sfc(1:3,max_blocks*npe))

    allocate(igrid_to_sfc(max_blocks))

    sfc=0
    allocate(Morton_start(0:npe-1),Morton_stop(0:npe-1))
    allocate(Morton_sub_start(0:npe-1),Morton_sub_stop(0:npe-1))

    allocate(nleafs_level(1:nlevelshi))

    allocate(coarsen(max_blocks,0:npe-1),refine(max_blocks,0:npe-1))
    coarsen=.false.
    refine=.false.
    if (nbufferx1/=0.or.nbufferx2/=0.or.nbufferx3/=0) then
       allocate(buffer(max_blocks,0:npe-1))
       buffer=.false.
    end if
    allocate(igrid_inuse(max_blocks,0:npe-1))
    igrid_inuse=.false.
    

    allocate(tree_root(1:ng1(1),1:ng2(1),1:ng3(1)))
    do ig3=1,ng3(1)
    do ig2=1,ng2(1)
    do ig1=1,ng1(1)
    nullify(tree_root(ig1,ig2,ig3)%node)
    end do
    end do
    end do

    ! define index ranges for ghost-cell swap
    call init_bc()

  end subroutine initialize_vars

  !> Set the index range over which the cell metrics are known.  With an even
  !> number of ghost cells this is just the ghosted block; with an odd number
  !> (or a staggered grid) ghost-cell prolongation needs one extra layer.
  subroutine set_geometry_index_ranges()
    use mod_global_parameters

    integer :: icase

    icase = mod(nghostcells,2)
    if (stagger_grid) icase = 1
    select case (icase)
    case (0)
       ixGextmin1=ixGlo1;ixGextmin2=ixGlo2;ixGextmin3=ixGlo3
       ixGextmax1=ixGhi1;ixGextmax2=ixGhi2;ixGextmax3=ixGhi3;
    case (1)
       ixGextmin1=ixGlo1-1;ixGextmin2=ixGlo2-1;ixGextmin3=ixGlo3-1
       ixGextmax1=ixGhi1+1;ixGextmax2=ixGhi2+1;ixGextmax3=ixGhi3+1;
    case default
       call mpistop("no such case")
    end select

  end subroutine set_geometry_index_ranges

  !> Allocate the cell metrics of all max_blocks blocks in one go.  The
  !> ixGext range is the one over which dvolume and ds are known, the face
  !> areas start one cell lower than ixG in every direction.
  subroutine alloc_geometry(geo, ixGmin1,ixGmin2,ixGmin3, ixGmax1,ixGmax2,&
       ixGmax3, ixGextmin1,ixGextmin2,ixGextmin3, ixGextmax1,ixGextmax2,&
       ixGextmax3)
    use mod_physicaldata, only: geo_t
    use mod_global_parameters, only: max_blocks, ndim

    type(geo_t), intent(inout) :: geo
    integer, intent(in)        :: ixGmin1,ixGmin2,ixGmin3,ixGmax1,ixGmax2,&
       ixGmax3, ixGextmin1,ixGextmin2,ixGextmin3,ixGextmax1,ixGextmax2,&
       ixGextmax3

    allocate( geo%x(ixGmin1:ixGmax1,ixGmin2:ixGmax2,ixGmin3:ixGmax3, 1:ndim,&
         1:max_blocks) )

    ! Everything past the positions is left unallocated in a Cartesian build.
    ! Nothing there reads it: the kernels that would (mod_finite_volume's flux
    ! divergence, mod_dt's cell sizes) are fypp-specialized and take their
    ! values from rnode instead, the AMR paths that would (prolong_2nd,
    ! coarsen_grid) sit in .not.slab_uniform branches that cannot run when
    ! GEOM is Cartesian, and get_volume_average has a slab branch.  The
    ! components stay declared in geo_t so that all of it still compiles; the
    ! matching state views are left unassociated by point_at_geometry, so a
    ! host reader that slips through faults instead of reading nonsense.
    !
    ! It is worth this much care because these are block-sized arrays for all
    ! max_blocks blocks a rank could ever own, allocated whether or not the
    ! blocks exist: on a 12 GB card, max_blocks = 4096 with block_nx = 16 is
    ! the difference between ~0.9 GB and ~4.1 GB of them.

  end subroutine alloc_geometry


end module mod_initialize