initialize_vars Subroutine

private subroutine initialize_vars()

Uses

  • proc~~initialize_vars~~UsesGraph proc~initialize_vars initialize_vars module~mod_amr_fct mod_amr_fct proc~initialize_vars->module~mod_amr_fct module~mod_fix_conserve mod_fix_conserve proc~initialize_vars->module~mod_fix_conserve module~mod_forest mod_forest proc~initialize_vars->module~mod_forest module~mod_geometry mod_geometry proc~initialize_vars->module~mod_geometry module~mod_ghostcells_update mod_ghostcells_update proc~initialize_vars->module~mod_ghostcells_update module~mod_global_parameters mod_global_parameters proc~initialize_vars->module~mod_global_parameters module~mod_mpi_wrapper mod_mpi_wrapper module~mod_amr_fct->module~mod_mpi_wrapper module~mod_fix_conserve->module~mod_mpi_wrapper module~mod_comm_lib mod_comm_lib module~mod_geometry->module~mod_comm_lib module~mod_ghostcells_update->module~mod_mpi_wrapper module~mod_basic_types mod_basic_types module~mod_global_parameters->module~mod_basic_types module~mod_connectivity mod_connectivity module~mod_global_parameters->module~mod_connectivity module~mod_constants mod_constants module~mod_global_parameters->module~mod_constants module~mod_physicaldata mod_physicaldata module~mod_global_parameters->module~mod_physicaldata module~mod_variables mod_variables module~mod_global_parameters->module~mod_variables mpi mpi module~mod_global_parameters->mpi module~mod_mpi_wrapper->mpi module~mod_variables->module~mod_basic_types

Initialize (and allocate) simulation and grid variables

Arguments

None

Calls

proc~~initialize_vars~~CallsGraph proc~initialize_vars initialize_vars proc~alloc_geometry alloc_geometry proc~initialize_vars->proc~alloc_geometry proc~check_pole_setup check_pole_setup proc~initialize_vars->proc~check_pole_setup proc~init_bc init_bc proc~initialize_vars->proc~init_bc proc~mpistop mpistop proc~initialize_vars->proc~mpistop proc~set_geometry_index_ranges set_geometry_index_ranges proc~initialize_vars->proc~set_geometry_index_ranges proc~set_pole set_pole proc~initialize_vars->proc~set_pole proc~check_pole_setup->proc~mpistop proc~init_bc->proc~mpistop mpi_abort mpi_abort proc~mpistop->mpi_abort proc~set_geometry_index_ranges->proc~mpistop proc~set_pole->proc~mpistop

Called by

proc~~initialize_vars~~CalledByGraph proc~initialize_vars initialize_vars proc~initialize_amrvac initialize_amrvac proc~initialize_amrvac->proc~initialize_vars proc~main main proc~main->proc~initialize_amrvac program~agile agile program~agile->proc~main

Source Code

  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