fill_geometry_device Subroutine

private subroutine fill_geometry_device(igrid)

Uses

  • proc~~fill_geometry_device~4~~UsesGraph proc~fill_geometry_device~4 fill_geometry_device module~mod_global_parameters mod_global_parameters proc~fill_geometry_device~4->module~mod_global_parameters 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_variables->module~mod_basic_types

Build block igrid's cell metrics directly on the device.

A uniform block is an analytic function of nothing but its corner and its spacing, both of which alloc_node has just put in rnode, so there is no reason to compute a dozen block-sized arrays on the CPU and ship them across the bus every time a node is (re)allocated - which, with AMR, is most regrids. Everything is built here instead, for whichever coordinate system the build selected, and nothing is handed back to the host - see sync_positions_host and sync_geometry_host in mod_geometry for that.

Two things about this routine are worth knowing before editing it.

First, the cell faces are primary and everything else is derived from them, via the RADIAL_CELL macro above. That is what lets one body serve both the uniformly spaced coordinate systems and the logarithmically stretched ones: a log-radial grid is uniform in xi = ln(r), so it is still exactly what this routine assumes - an analytic function of the block's corner and a constant spacing - and only the step from logical faces to physical ones differs.

Second, the position written into bgeo%x is the volume barycentre of the cell, not the midpoint of its faces. The solution is stored as a cell average, and a cell average equals the point value at the barycentre to second order, because the linear term of the Taylor expansion integrates to zero only about that point. So the barycentre is where the state, the initial condition and any position-dependent source term belong. The consequence to remember is that x +/- dx/2 is not a cell face any more: the two places that need a face position rebuild it from rnode instead (the face coordinates in mod_finite_volume, and calc_x's corner grid). Only the radial direction, and spherical's polar direction, are affected - phi and cylindrical z carry uniform weight, so their barycentres are their midpoints.

x uses the two-sided form (measured from rnode's lower corner in the mesh and its upper corner in the outer ghost layer, so the overlapping ghost cells of neighbouring blocks agree to the last bit), while the volumes and cell sizes use the one-sided form over the extended range that the old code called xext.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: igrid

Called by

proc~~fill_geometry_device~4~~CalledByGraph proc~fill_geometry_device~4 fill_geometry_device proc~alloc_node~4 alloc_node proc~alloc_node~4->proc~fill_geometry_device~4

Source Code

  subroutine fill_geometry_device(igrid)
    use mod_global_parameters

    integer, intent(in) :: igrid

    integer          :: ix1,ix2,ix3
    integer          :: ixCoGmax1,ixCoGmax2,ixCoGmax3
    ! last cell in each direction whose centre is measured from the lower corner
    integer          :: iup1,iup2,iup3
    ! cell spacing of the block and of its coarse representative
    double precision :: d1,d2,d3, c1,c2,c3
    double precision :: xmin1,xmin2,xmin3, xmax1,xmax2,xmax3
    ! logical cell centre (p) and extended-range logical cell centre (e)
    double precision :: p1,p2,p3, e1,e2

    d1=rnode(rpdx1_,igrid); d2=rnode(rpdx2_,igrid); d3=rnode(rpdx3_,igrid)
    c1=2.0d0*d1; c2=2.0d0*d2; c3=2.0d0*d3
    xmin1=rnode(rpxmin1_,igrid); xmin2=rnode(rpxmin2_,igrid)
    xmin3=rnode(rpxmin3_,igrid)
    xmax1=rnode(rpxmax1_,igrid); xmax2=rnode(rpxmax2_,igrid)
    xmax3=rnode(rpxmax3_,igrid)

    iup1=ixMhi1-nghostcells; iup2=ixMhi2-nghostcells; iup3=ixMhi3-nghostcells

    ixCoGmax1=(ixGhi1-2*nghostcells)/2+2*nghostcells
    ixCoGmax2=(ixGhi2-2*nghostcells)/2+2*nghostcells
    ixCoGmax3=(ixGhi3-2*nghostcells)/2+2*nghostcells

    ! ---- positions, and (curvilinear only) the face areas built from them ----
    ! The three "if (ix==ixGlo)" branches fill the lower face of the first
    ! cell in each direction, which the flux update reaches as surfaceC(ix-1,.)
    ! That lower face is just fL of the first cell, which the loop already has.
    ! A Cartesian build has no face areas to fill: its flux update divides by
    ! the uniform rnode spacing instead, so surfaceC is never allocated.
    ! The cell centre is measured from the block's lower corner inside the
    ! mesh, and from its upper corner in the outer ghost layer, so that the
    ! overlapping ghost cells of two neighbouring blocks come out identical to
    ! the last bit.
     
    do ix3=ixGlo3,ixGhi3
       do ix2=ixGlo2,ixGhi2
          do ix1=ixGlo1,ixGhi1
             if (ix1 <= iup1) then
                p1 = xmin1 + (dble(ix1-nghostcells)-half)*d1
             else
                p1 = xmax1 + (dble(ix1-ixMhi1)-half)*d1
             end if
             if (ix2 <= iup2) then
                p2 = xmin2 + (dble(ix2-nghostcells)-half)*d2
             else
                p2 = xmax2 + (dble(ix2-ixMhi2)-half)*d2
             end if
             if (ix3 <= iup3) then
                p3 = xmin3 + (dble(ix3-nghostcells)-half)*d3
             else
                p3 = xmax3 + (dble(ix3-ixMhi3)-half)*d3
             end if

             ! uniform weight in every direction, so the barycentre of a cell
             ! is its midpoint and the logical coordinate is the physical one
             bgeo%x(ix1,ix2,ix3,1,igrid) = p1
             bgeo%x(ix1,ix2,ix3,2,igrid) = p2
             bgeo%x(ix1,ix2,ix3,3,igrid) = p3
          end do
       end do
    end do


    ! ---- the one-level-coarser representative ----
    ! Its positions are measured from the lower corner throughout - there is no
    ! neighbouring coarse block whose ghost cells have to match - and its
    ! spacing is doubled, so one loop covers positions, volumes and areas.
    ! As above, only the positions exist in a Cartesian build.
     
    do ix3=1,ixCoGmax3
       do ix2=1,ixCoGmax2
          do ix1=1,ixCoGmax1
             p1 = xmin1 + (dble(ix1-nghostcells)-half)*c1
             p2 = xmin2 + (dble(ix2-nghostcells)-half)*c2
             p3 = xmin3 + (dble(ix3-nghostcells)-half)*c3

             bgeoc%x(ix1,ix2,ix3,1,igrid) = p1
             bgeoc%x(ix1,ix2,ix3,2,igrid) = p2
             bgeoc%x(ix1,ix2,ix3,3,igrid) = p3
          end do
       end do
    end do

    ! Nothing is copied back here. alloc_node runs on every regrid, mostly for
    ! blocks no host routine will look at, so the fetch belongs at the sites
    ! that actually read the geometry: sync_positions_host before host code
    ! that needs ps(igrid)%x, sync_geometry_host before output.

  end subroutine fill_geometry_device