setaux default routines

For GeoClaw, see setaux routine in GeoClaw.

Below are the default setaux library routines for Classic and AMRClaw. By default these do nothing. If you wish to specify aux arrays you will need to copy one of these files to your application directory and modify it as needed. Remember to change to Makefile to point to the proper version.

$CLAW/classic/src/1d/setaux.f90:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
subroutine setaux(mbc,mx,xlower,dx,maux,aux)

    ! Called at start of computation before calling qinit, and
    ! when AMR is used, also called every time a new grid patch is created.
    ! Use to set auxiliary arrays aux(1:maux, 1-mbc:mx+mbc, 1-mbc:my+mbc).
    ! Note that ghost cell values may need to be set if the aux arrays
    ! are used by the Riemann solver(s).
    !
    ! This default version does nothing. 
 
    implicit none
    integer, intent(in) :: mbc,mx,maux
    real(kind=8), intent(in) :: xlower,dx
    real(kind=8), intent(out) ::  aux(maux,1-mbc:mx+mbc)

end subroutine setaux

$CLAW/classic/src/2d/setaux.f90:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
subroutine setaux(mbc,mx,my,xlower,ylower,dx,dy,maux,aux)

    ! Called at start of computation before calling qinit, and
    ! when AMR is used, also called every time a new grid patch is created.
    ! Use to set auxiliary arrays aux(1:maux, 1-mbc:mx+mbc, 1-mbc:my+mbc).
    ! Note that ghost cell values may need to be set if the aux arrays
    ! are used by the Riemann solver(s).
    !
    ! This default version does nothing. 
 
    implicit none
    integer, intent(in) :: mbc,mx,my,maux
    real(kind=8), intent(in) :: xlower,ylower,dx,dy
    real(kind=8), intent(out) ::  aux(maux,1-mbc:mx+mbc,1-mbc:my+mbc)

end subroutine setaux

$CLAW/classic/src/3d/setaux.f90:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
subroutine setaux(mbc,mx,my,mz,xlower,ylower,zlower,dx,dy,dz,maux,aux)

    ! Called at start of computation before calling qinit, and
    ! when AMR is used, also called every time a new grid patch is created.
    ! Use to set auxiliary arrays 
    !   aux(1:maux, 1-mbc:mx+mbc, 1-mbc:my+mbc, 1-mbc:mz+mbc).
    ! Note that ghost cell values may need to be set if the aux arrays
    ! are used by the Riemann solver(s).
    !
    ! This default version does nothing. 
 
    implicit none
    integer, intent(in) :: mbc,mx,my,mz,maux
    real(kind=8), intent(in) :: xlower,ylower,zlower,dx,dy,dz
    real(kind=8), intent(out) ::  aux(maux,1-mbc:mx+mbc,1-mbc:my+mbc)

end subroutine setaux

setaux routine in GeoClaw

In GeoClaw, there is a library routine that sets aux(1,i,j) to the cell average of the bathymetry, aux(2,i,j) to the ratio of the true cell area to dx*dy (the capacity function), and aux(3,i,j) to the length ratio of the bottom edge to dx (The latter two quantities vary with latitude when coordinate_system == 2).

$CLAW/geoclaw/src/2d/shallow/setaux.f90:

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
! Set auxiliary arrays
!
! Routine that sets the values found in the aux array including topography and
! other fields that are used in the source terms and other things.
!
! In the default routine sets the following fields depending on the input
! parameters present:
!  (1) topography
!  (2) capacity (set if coordinate_system == 2)
!  (3) length ratio of edge (set if coordinate_system == 2)
!  (frition_index). Location of manning's N (variable_friction)
!  (wind_index, wind_index + 1).  Location of x and y wind speeds (wind_forcing)
!  (pressure_index).  Location of pressure field (pressure_forcing)
!
! There are a couple of additional things of note:
!
!  - This routine is called anytime a grid is created during re-gridding to fill
!    in the aux arrays.
!  - Built-in to this routine is an ability to copy aux values from grids from
!    the same level.  If you do not want this behavior you will need to modify
!    this routine.
!  - If you are using periodic BCs then you must handle this here.  Look to the
!    setting of topography for an example of how this may need to be done.
!  - If your aux arrays are time dependent and set somewhere else it may be
!    prudent to set a default value here as is done with the wind and pressure.
!
subroutine setaux(mbc,mx,my,xlow,ylow,dx,dy,maux,aux)

    use amr_module, only: mcapa, xupper, yupper, xlower, ylower, NEEDS_TO_BE_SET
    use amr_module, only: xperdom, yperdom

    use geoclaw_module, only: coordinate_system, earth_radius, deg2rad
    use geoclaw_module, only: sea_level, ambient_pressure

    use storm_module, only: wind_forcing, pressure_forcing
    use storm_module, only: wind_index, pressure_index, set_storm_fields

    use friction_module, only: variable_friction, friction_index
    use friction_module, only: set_friction_field

    use topo_module

    use adjoint_module, only : adjoint_flagging,innerprod_index

    implicit none

    ! Arguments
    integer, intent(in) :: mbc,mx,my,maux
    real(kind=8), intent(in) :: xlow,ylow,dx,dy
    real(kind=8), intent(inout) :: aux(maux,1-mbc:mx+mbc,1-mbc:my+mbc)

    ! Locals
    integer :: ii,jj,m, iint,jint
    real(kind=8) :: x,y,xm,ym,xp,yp,topo_integral
    real(kind=8) :: xper,yper,xperm,yperm,xperp,yperp
    character(len=*), parameter :: aux_format = "(2i4,4d15.3)"
    integer :: skipcount,iaux,ilo,jlo

    ! Lat-Long coordinate system in use, check input variables
    if (coordinate_system == 2) then
        if (mcapa /= 2 .or. maux < 3) then
            print *,'ERROR in setaux:  for coordinate_system==2'
            print *,'     need mcapa == 2 and maux >= 3'
            print *,'     have mcapa = ',mcapa,'  maux = ',maux
            stop
        endif
    endif

    ! Compute integer indices based off same corner of domain to reduce round
    ! off discrepancies
    ilo = floor((xlow - xlower + .05d0*dx)/dx)
    jlo = floor((ylow - ylower + .05d0*dy)/dy)

    ! Set geometry values
    if (coordinate_system == 1) then
        if (maux == 0 .and. mcapa > 0) then
            print *, "ERROR:  Capacity array requested but number of aux"
            print *, "variables is set to 0."
            stop
        end if
    else if (coordinate_system == 2) then
        do jj = 1 - mbc, my + mbc
            do ii = 1 - mbc, mx + mbc
                ym = ylower + (jlo+jj-1.d0) * dy
                yp = ylower + (jlo+jj) * dy
                aux(2,ii,jj) = deg2rad * earth_radius**2                      &
                                * (sin(yp * deg2rad) - sin(ym * deg2rad)) / dy
                aux(3,ii,jj) = ym * deg2rad
            end do
        end do
    end if

    ! Storm fields if used
    if (wind_forcing) then
        aux(wind_index, :, :) = 0.d0
        aux(wind_index + 1, :, :) = 0.d0
    endif
    if (pressure_forcing) then
        aux(pressure_index, :, :) = ambient_pressure
    endif

    ! Innerproduct field if used
    if (adjoint_flagging) then
        do jj=1-mbc,my+mbc
            do ii=1-mbc,mx+mbc
                if (aux(1,ii,jj) .eq. NEEDS_TO_BE_SET) then
                    aux(innerprod_index,ii,jj) = 0.d0
                endif
            enddo
        enddo
    endif


    ! ==========================================================================
    !  Set Bathymetry
    ! Set analytical bathymetry here if requested
    if (test_topography > 0) then
        do jj = 1 - mbc, my + mbc
            ym = ylower + (jlo+jj-1.d0) * dy
            yp = ylower + (jlo+jj) * dy
            y = 0.5d0*(ym+yp)
            do ii = 1 - mbc, mx + mbc
                xm = xlower + (ilo+ii-1.d0) * dx
                xp = xlower + (ilo+ii) * dx
                x = 0.5d0*(xm+xp)
                aux(1,ii,jj) = test_topo(x, y)
            end do
        end do

    ! Set bathymetry based on input files
    else if (mtopofiles > 0) then

        do jj=1-mbc,my+mbc
            ym = ylower + (jlo+jj-1.d0) * dy
            yp = ylower + (jlo+jj) * dy
            y = 0.5d0*(ym+yp)


            do ii=1-mbc,mx+mbc
                xm = xlower + (ilo+ii-1.d0) * dx
                xp = xlower + (ilo+ii) * dx
                x = 0.5d0*(xm+xp)

                ! Parameter NEEDS_TO_BE_SET initialized in amr_module.f90
                ! saves time by otherwise copying instead of reinitializing
                if (aux(1,ii,jj) .ne. NEEDS_TO_BE_SET) then
                    cycle
                endif

                topo_integral = 0.d0
                if ((y>yupper).or.(y<ylower).or.(x>xupper).or.(x<xlower)) then
                    if (.not.(xperdom .or. yperdom)) then
                        ! Skip setting as this cell sticks out of the physical
                        ! domain and we are not setting periodic BCs
                        cycle

                    else
                        ! We evaluate the periodic BC topography by computing
                        ! the appropriate periodic grid cell coordinates and
                        ! again calling cellgridintegrate
                        call wrap_coords(x,y,xperm,xper,xperp,yperm,yper, &
                                         yperp,dx,dy)
                        call cellgridintegrate(topo_integral,  &
                            xperm,xper,xperp,yperm,yper,yperp, &
                            xlowtopo,ylowtopo,xhitopo,yhitopo,dxtopo,dytopo, &
                            mxtopo,mytopo,mtopo,i0topo,mtopoorder, &
                            mtopofiles,mtoposize,topowork)
                    endif
                else
                    ! Cell does not extend outside of physical domain
                    call cellgridintegrate(topo_integral,xm,x,xp,ym,y,yp, &
                            xlowtopo,ylowtopo,xhitopo,yhitopo,dxtopo,dytopo, &
                            mxtopo,mytopo,mtopo,i0topo,mtopoorder, &
                            mtopofiles,mtoposize,topowork)
                endif

                ! Correct for geometry
                if (coordinate_system == 2) then
                    aux(1,ii,jj) = topo_integral / (dx * dy * aux(2,ii,jj))
                else
                    aux(1,ii,jj) = topo_integral / (dx * dy)
                endif
            enddo
        enddo
    else
        print *, "ERROR:  There is no way to set bathymetry!  Either "
        print *, "        provide topography files or request topography "
        print*,  "        defined by a function."
        stop
    end if

    ! Copy topo to ghost cells if outside physical domain and not periodic
    if (.not. yperdom) then
        do jj=1-mbc,my+mbc
            y = ylower + (jlo+jj-.5d0) * dy
            if ((y < ylower) .or. (y>yupper)) then
                do ii=1-mbc,mx+mbc
                    x = xlower + (ilo+ii-.5d0) * dx
                    iint = ii + max(0, ceiling((xlower-x)/dx)) &
                         - max(0, ceiling((x-xupper)/dx))
                    jint = jj + max(0, ceiling((ylower-y)/dy)) &
                         - max(0, ceiling((y-yupper)/dy))
                    aux(1,ii,jj) = aux(1,iint,jint)
                enddo
            endif
        enddo
    endif
    if (.not. xperdom) then
        do ii=1-mbc,mx+mbc
            x =  xlower + (ilo+ii-.5d0) * dx
            if ((x < xlower) .or. (x > xupper)) then
                do jj=1-mbc,my+mbc
                    y = ylower + (jlo+jj-.5d0) * dy
                    iint = ii + max(0, ceiling((xlower-x)/dx)) &
                         - max(0, ceiling((x-xupper)/dx))
                    jint = jj + max(0, ceiling((ylower-y)/dy)) &
                         - max(0, ceiling((y-yupper)/dy))
                    aux(1,ii,jj) = aux(1,iint,jint)
                enddo
            endif
        enddo
    endif

    ! If using a variable friction field initialize the coefficients to 0
    if (variable_friction) then
        call set_friction_field(mx, my, mbc, maux, xlow, ylow, dx, dy, aux)
    endif

contains

    ! Provide wrapper function for providing periodic coordinates
    subroutine wrap_coords(x,y,xperm,xper,xperp,yperm,yper,yperp,dx,dy)

        use amr_module, only: xperdom, yperdom, xupper, yupper, xlower, ylower

        implicit none

        ! Arguments
        real(kind=8), intent(in) :: x,y,dx,dy
        real(kind=8), intent(out) :: xperm,xper,xperp,yperm,yper,yperp
        real(kind=8) :: xdomain, ydomain

        ! need size of domain for wrapping
        xdomain = xupper - xlower
        ydomain = yupper - ylower
        xper = x
        yper = y

        ! test for wrapping of coordinates
        if (x > xupper .and. xperdom) xper = xper - xdomain
        if (x < xlower .and. xperdom) xper = xper + xdomain
        if (y > yupper .and. yperdom) yper = yper - ydomain
        if (y < ylower .and. yperdom) yper = yper + ydomain

        ! adjust rest of variables
        xperm = xper - 0.5d0 * dx
        xperp = xper + 0.5d0 * dx
        yperm = yper - 0.5d0 * dy
        yperp = yper + 0.5d0 * dy

    end subroutine wrap_coords

end subroutine setaux