In this tutorial you will be led through the process of developing a simple 2D diffusion solver using the PARAMESH package. Follow the steps outlined below. When you are finished you will have created and run a simple calculation with an adaptive mesh.
As explained in the introduction, PARAMESH enables you to build your application with or without permanently allocated memory for storing guardcell data. This first tutorial assumes you will allocate permanent guardcell storage space. A second tutorial which can be accessed from the bottom of this page takes you through the steps required to build the same application without permanent guardcell storage. A third tutorial is also available from the bottom of this page which takes you through the process of using PARAMESH as a library. Finally, a fourth tutorial is available on the following pages which describes how to add flux conservation to the previous tutorials.
The simple problem we tackle here is to integrate a 2D diffusion
equation.
We will use grid blocks of size 4x4 with 1 guard cell at each boundary.
The data points defining U are assumed to be known at cell centers.
We use a very simple explicit time integration scheme, and a 4-pt
centered
second order accurate difference to represent the 
operator, ie.
U(i,j,t+dt) = U(i,j,t) + dt * A /(dx*dx)
A = U(i+1,j,t) + U(i-1,j,t) + U(i,j+1,t) + U(i,j-1,t) - 4*U(i,j,t)
In the tutorial described on this page, we will assume that you will be using MPI. You must have a version of MPI installed on your system to use this tutorial. You can get a freely available version of MPI which can be installed on most systems from the MPICH web site
Step 1. - Preparations
Step 2. - Edit the header file paramesh_preprocessor.fh.
#define REAL8
Note, this tells paramesh that you will be using real*8 as the default
for reals, but it does NOT pass this information to the compiler !
If you choose double precision, you must modify the makefile so
that the compiler also knows !
!#define VAR_DT
!#define PRED_CORR
!#define EMPTY_CELLS
#define DIAGONALS
#define N_DIM 2
!#define NO_PERMANENT_GUARDCELLS
!#define ADVANCE_ALL_LEVELS
#define NX_B 4
#define NY_B 4
#define MAX_BLOCKS 100
#define N_GUARD_CELLS 1
#define N_GUARD_CELLS_WORK 1
#define N_VAR 1
#define N_FACEVAR 0
#define N_VAR_EDGE 0
#define N_VAR_CORN 0
#define N_VAR_WORK 1
#define N_FLUX_VAR 1
#define N_EDGE_VAR 0
All other settings in this file can be left with their default values.
Step 3. - Create and Edit the tutorial's main makefile.
Explanation:
This makefile will execute makefiles in each of the paramesh source
sub-directories (/headers, /source, and if using mpi in /mpi_source).
Step 4. - Edit the tutorial sub-makefile.
main := tutorial.F90
sources := amr_1blk_bcset.F90
CMD = tutor
Explanation:
To begin with we will work on just 1 processor. The file tutorial.F90
now
contains a template main program,
and the file amr_1blk_bcset.F90 a template of the boundary condition
routine, which we need to provide to PARAMESH. We will edit these
templates
soon.
The executable which will be made will be called 'tutor'.
Step 5. - Modify the main program template.
Explanation:
Section 1 is a list of include files and definitions of some local
variables.
Section 2 initializes the amr package, and reports the number of active
processors. Section 3 establishes the initial grid. It does this in the
simplest possible way. First 1 grid block is established on processor 0
which covers the entire computational domain. This is marked for
refinement. Then we enter a loop which we repeat twice. Inside this
loop we mark all existing blocks for refinement, and then implement
these refinements. The first pass through refines the first block, and
creates its 4 children. On the second pass
through, these 4 children and marked for refinement, and then refined.
The result will be a tree of grid blocks with 16 leaf blocks at
refinement level 3, 4 parents at refinement level 2, and the original
block at level 1.
Print statements have been provided to output details of the grid
blocks
which have been set up.
Step 6. - Modify the file amr_1blk_bcset.F90
! if(ibc.eq. ???? ) thenand its corresponding endif . Change the ???? in the if statement to any integer less than or equal to -20.
! unk1(:,i,j,k,idest) = ???? !<<<<< USER EDITand replace the right hand side of this line with 0.0
Explanation:
This routine is needed as part of the mechanism by which the User
supplies any non-periodic boundary conditions which might be required.
In our example, we are using periodic boundary conditions. These are
implemented
by identifying blocks at the left boundary as neighbors of blocks at
the right
boundary, and vice versa. Because we have no blocks whose neighbors
addresses will satisfy
the ibc if test in this example, this routine will do nothing.
We include this step in the tutorial merely to introduce you to the
mechanism you must use when setting up aperiodic boundary conditions.
Step 7. - Build the executable.
gmake -f make_tutor your_tutorial
Step 8. - Run the executable.
mpirun -np 1 your_tutorial/tutor
Congratulations, you have now begun to do amr !
If everything went according to plan you should have generated a short output listing which concludes with something equivalent to the following lines (the order in which the blocks are listed may vary slightly, from one machine to another):
pe / blk / blk-coords / blk-sizes 0, 1, 2*0.E+0, 2*8. 0, 2, 2*-3., 2*2. 0, 3, 2*-2., 2*4. 0, 4, -1., -3., 2*2. 0, 5, -3., -1., 2*2. 0, 6, 2*-1., 2*2. 0, 7, 2., -2., 2*4. 0, 8, 1., -3., 2*2. 0, 9, 3., -3., 2*2. 0, 10, 1., -1., 2*2. 0, 11, 3., -1., 2*2. 0, 12, -3., 1., 2*2. 0, 13, -2., 2., 2*4. 0, 14, -1., 1., 2*2. 0, 15, -3., 3., 2*2. 0, 16, -1., 3., 2*2. 0, 17, 2*2., 2*4. 0, 18, 2*1., 2*2. 0, 19, 3., 1., 2*2. 0, 20, 1., 3., 2*2. 0, 21, 2*3., 2*2.
The full grid is shown here with block boundaries drawn in black and grid cells outlined in red.
This is a list of all the grid blocks that you have created. Since we are currently only using 1 processor, all the blocks are listed as being on processor 0. The second integer on each line is the block address on processor 0, the next 2 floating point numbers specify the x and y coordinates of the center of each block, and the last two numbers are the block lengths along the x and y axes.
Step 9. - Initializing the solution.
Step 10. - Edit amr_initial_soln.F90
dx = bsize(1,lb)/real(nxb)
dy = bsize(2,lb)/real(nyb)
unk(1,i,j,k,lb) = 1.0
xi = bnd_box(1,1,lb) + dx*(real(i-nguard0)-.5)
yi = bnd_box(1,2,lb) + dy*(real(j-nguard0)-.5)
if( abs(xi).lt.1.0 .and. abs(yi).lt.1.0) then
unk(1,i,j,k,lb) = 10.0
endif
Explanation:
This step has modified the routine amr_initial_soln to generate the
initial solution we want.
Step 11. - Edit tutorial.F
do lb=1,lnblocks
if(coord(1,lb).eq.1.0.and.coord(2,lb).eq.1.0) then
do j=1,nyb+2*nguard
write(*,50) j,(unk(1,i,j,1,lb),i=1,nxb+2*nguard)
50 format(1x,i3,6(2x,f7.4))
enddo
endif
enddo
Explanation:Step 12. - Remake and run.
gmake -f make_tutor your_tutorial
mpirun -np 1 your_tutorial/tutor
You have now initialized the solution array unk(1,:,:,:,:) on all the grid blocks of the initial grid. As proof, the last six lines of your output show the data values on the centered at (1.0,1.0). It should look like this:
1 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 2 0.0000 10.0000 10.0000 1.0000 1.0000 0.0000 3 0.0000 10.0000 10.0000 1.0000 1.0000 0.0000 4 0.0000 1.0000 1.0000 1.0000 1.0000 0.0000 5 0.0000 1.0000 1.0000 1.0000 1.0000 0.0000 6 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
This block is located at 0 < x < 2 and 0 < y < 2. It straddles one corner of the high density region. Notice, the 4x4 block interior has been initialized with non-zero values and there is a layer of guard cells surrounding the block which are currently all set to 0.0.
The complete initial state is shown here, with the block boundaries
superimposed
in black and the grid cells outlined in red.
Step 13. - Filling Guardcells
Edit tutorial.F90
do lb=1,lnblocks
if(coord(1,lb).eq.1.0.and.coord(2,lb).eq.1.0) then
do j=1,nyb+2*nguard
write(*,50) j,(unk(1,i,j,1,lb),i=1,nxb+2*nguard)
50 format(1x,i3,6(2x,f7.4))
enddo
endif
enddo
Step 14. - Remake and run.
gmake -f make_tutor your_tutorial
mpirun -np 1 your_tutorial/tutor
Now the last 6 lines which report the solution array in our selected block looks like this.
1 10.0000 10.0000 10.0000 1.0000 1.0000 1.0000 2 10.0000 10.0000 10.0000 1.0000 1.0000 1.0000 3 10.0000 10.0000 10.0000 1.0000 1.0000 1.0000 4 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 5 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 6 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000Notice, the guard cell layer has been filled with the correct data from the neighboring blocks.
Step 15. - Constructing a routine to test refinement levels.
Explanation:
This routines function is to set the logical arrays 'refine' and
'derefine'
to be true or false for each grid block. It does this by computing a
local
error measure and testing to see if the acceptable limits on this error
measure are violated anywhere within each leaf block. There is one
subtle
feature of this example which we should mention. The data used as input
for
the error calculation is first placed in the 'work' array. This is a
special cell centered data-structure, because, unlike 'unk', when the
guardcell filling routines operate on 'work' they limit their scope to
1 data word per grid cell. 'work' is dimensioned inside the header file
workspace.fh. It can have a number of guard cell layers set differently
from that of the solution data-structure. By using 'work' here, and
setting
its guard cell number to be larger than the number of layers
guarding the solution data-structure, we can help grid blocks to refine
in advance of the arrival of disturbances.
We will use a very simple error measure here, based on variation of the solution between neighboring grid cells.
error(:,:,:) = 0.
do k=klw,kuw
do j=jlw+1,juw-1
do i=ilw+1,iuw-1
error1 = abs(work(i+1,j,k,lb,1)-work(i,j,k,lb,1))
error2 = abs(work(i-1,j,k,lb,1)-work(i,j,k,lb,1))
error3 = abs(work(i,j+1,k,lb,1)-work(i,j,k,lb,1))
error4 = abs(work(i,j-1,k,lb,1)-work(i,j,k,lb,1))
error_num = max( error1,error2,error3,error4 )
error_den = max( work(i,j,k,lb,1) ,work(i+1,j,k,lb,1), &
work(i-1,j,k,lb,1),work(i,j+1,k,lb,1), &
work(i,j-1,k,lb,1), 1.0e-6 )
error(i,j,k) = error_num/error_den
enddo
enddo
enddo
Step 16. - Edit tutorial.F90
lrefine_max = 4
if(mype.eq.0) write(*,*) ' pe blk refine derefine', &
' curr.ref.level'
call MPI_BARRIER(MPI_COMM_WORLD, ierr)
do l=1,lnblocks
write(*,51) mype,l,refine(l),derefine(l),lrefine(l)
51 format(1x,i3,2x,i3,2x,l8,2x,l8,10x,i3)
enddo
if(mype.eq.0) write(*,*) 'pe / blk / blk-coords / blk-sizes'
call MPI_BARRIER(MPI_COMM_WORLD, ierr)
do l=1,lnblocks
write(*,*) mype,l,(coord(i,l),i=1,ndim),(bsize(j,l),j=1,ndim)
enddo
Step 17. - Set the number of guard cell layers for 'work'
Step 18. - Remake and run.
gmake -f make_tutor your_tutorial
mpirun -np 1 your_tutorial/tutor
Explanation:
The changes that you have just made, analyzed the error estimate on
each
existing grid block, and marked some blocks for additional refinement.
In your new output there will be a section looking like this
(the order of lines may be slightly different) :
pe blk refine derefine curr.ref.level 0 1 F F 1 0 2 F T 3 0 3 T F 2 0 4 F T 3 0 5 F T 3 0 6 T F 3 0 7 T F 2 0 8 F T 3 0 9 F T 3 0 10 T F 3 0 11 F T 3 0 12 F T 3 0 13 T F 2 0 14 T F 3 0 15 F T 3 0 16 F T 3 0 17 T F 2 0 18 T F 3 0 19 F T 3 0 20 F T 3 0 21 F T 3This is telling you that the test in amr_test_refinement marked blocks 3, 6, 7, 10, 13, 14, 17, and 18 for further refinement. However blocks 3, 7, 13, and 17 are parent blocks at level 2 and so their refinement flags will be ignored. Blocks 6, 10, 14 and 18 will be refined. Notice also that blocks 2,4,5,8,9,11,12,15,16,19,20 and 21 have been marked for derefinement. However each of these blocks has a sibling which has not been marked for derefinement ( in fact all their siblings have been marked for refinement ), and so these particular derefinement choices will be cancelled by PARAMESH. You can see that this is the case in the listing of the new grid which follows. The next section of output reports some details of the refinement procedure, and the final section is a list of the new grid, which now has 37 blocks in total.
pe / blk / blk-coords / blk-sizes 0, 1, 2*0.E+0, 2*8. 0, 2, 2*-3., 2*2. 0, 3, 2*-2., 2*4. 0, 4, -1., -3., 2*2. 0, 5, -3., -1., 2*2. 0, 6, 2*-1.5, 2*1. 0, 7, 2*-1., 2*2. 0, 8, -0.5, -1.5, 2*1. 0, 9, -1.5, -0.5, 2*1. 0, 10, 2*-0.5, 2*1. 0, 11, 2., -2., 2*4. 0, 12, 1., -3., 2*2. 0, 13, 3., -3., 2*2. 0, 14, 0.5, -1.5, 2*1. 0, 15, 1., -1., 2*2. 0, 16, 1.5, -1.5, 2*1. 0, 17, 0.5, -0.5, 2*1. 0, 18, 1.5, -0.5, 2*1. 0, 19, 3., -1., 2*2. 0, 20, -3., 1., 2*2. 0, 21, -2., 2., 2*4. 0, 22, -1.5, 0.5, 2*1. 0, 23, -1., 1., 2*2. 0, 24, -0.5, 0.5, 2*1. 0, 25, -1.5, 1.5, 2*1. 0, 26, -0.5, 1.5, 2*1. 0, 27, -3., 3., 2*2. 0, 28, -1., 3., 2*2. 0, 29, 2*1., 2*2. 0, 30, 2*0.5, 2*1. 0, 31, 2*2., 2*4. 0, 32, 1.5, 0.5, 2*1. 0, 33, 0.5, 1.5, 2*1. 0, 34, 2*1.5, 2*1. 0, 35, 3., 1., 2*2. 0, 36, 1., 3., 2*2. 0, 37, 2*3., 2*2.The complete new grid is illustrated below.
Step 19. - Create the routine to update the solution.
subroutine advance_soln(mype,time,dt)
end subroutine advance_soln
integer :: mype
real :: time,dt
real old_soln(il_bnd:iu_bnd,jl_bnd:ju_bnd,kl_bnd:ku_bnd)
making sure that they appear after the use statements.
call amr_timestep(dt,dtmin,dtmax,mype)
old_soln(:,:,:) = unk(1,:,:,:,lb)
do k=kl_bnd+nguard*k3d,ku_bnd-nguard*k3d
do j=jl_bnd+nguard*k2d,ju_bnd-nguard*k2d
do i=il_bnd+nguard,iu_bnd-nguard
unk(1,i,j,k,lb) = old_soln(i,j,k) + dt/(dx*dx)* ( &
old_soln(i+1,j,k) + old_soln(i-1,j,k) + &
old_soln(i,j+1,k) + old_soln(i,j-1,k) - &
old_soln(i,j,k)*4.0 )
enddo
enddo
enddo
time = time + dt
Step 20. - Create the timestep routine.
rho => unk(1,:,:,:,l)
vx => unk(2,:,:,:,l)
vy => unk(3,:,:,:,l)
vz => unk(4,:,:,:,l)
real, parameter :: courant=.1, kappa=1.0
dtl = courant*dx*dx/kappa
Step 21. - Modify main program to call the advance_soln routine.
write(*,*) 'dt = ',dt
do lb=1,lnblocks
if(coord(1,lb).eq.1.0.and.coord(2,lb).eq.1.0) then
do j=1,nyb+2*nguard
write(*,50) j,(unk(1,i,j,1,lb),i=1,nxb+2*nguard)
enddo
endif
enddo
Step 22. - Remake and run.
gmake -f make_tutor your_tutorial
mpirun -np 1 your_tutorial/tutor
Explanation:
You have now advanced the solution through 1 timestep, and the output
section immediately after the call to advance_soln will show how the
data on the block centered on (1.0,1.0) has been diffused. The data
should look like this:
dt = 2.500000037E-2 1 10.0000 10.0000 10.0000 1.0000 1.0000 1.0000 2 10.0000 10.0000 9.1000 1.9000 1.0000 1.0000 3 10.0000 9.1000 8.2000 1.9000 1.0000 1.0000 4 1.0000 1.9000 1.9000 1.0000 1.0000 1.0000 5 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 6 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000Note, at this point the cell interior (indeces 2-5 in both x and y) are correct, but the guardcells (indeces 1 and 6) have not yet been updated.
After the solution has been advanced on the block interiors, we test the solution to see if refinement is required. In this case refinement is selected for the 4 blocks around the center of the domain. These are refined, and the solution is prolonged to the newly created blocks there. The complete updated solution after these steps is shown here.
Step 23. - Run for 250 timesteps.
if(mype.eq.0) write(*,*) 'iteration ',istep, &
' no of blocks = ',lnblocks
if(mype.eq.0) write(*,*) 'pe / blk / blk-coords / blk-sizes'
call MPI_BARRIER(MPI_COMM_WORLD, ierr)
do l=1,lnblocks
write(*,*) mype,l,(coord(i,l),i=1,ndim),(bsize(j,l),j=1,ndim)
enddo
gmake -f make_tutor your_tutorial
mpirun -np 1 your_tutorial/tutor
Explanation:
In your output you will notice that before the first timestep we had
21 blocks, with uniform refinement at level 3 throughout the
computational
domain. On the first timestep (iteration 1) the 4 blocks at the center
containing the high data values were refined, adding 16 blocks to make
a total of 37. After the seventeenth timestep the solution has diffused
outward
so that all the outer level 3 blocks, except for those on the corners,
are now all marked for refinement, adding another 32 child blocks at
level 4,
for a total of 69.
After 250 timesteps, we can see from the final block listing below that block number 53 is a leaf block located near the origin (it has coordinates x=0.5, y=0.5 ; the line order may be slightly different in your output).
pe / blk / blk-coords / blk-sizes 0, 1, 2*-3., 2*2. 0, 2, 2*-2., 2*4. 0, 3, 2*0.E+0, 2*8. 0, 4, -1., -3., 2*2. 0, 5, -1.5, -3.5, 2*1. 0, 6, -0.5, -3.5, 2*1. 0, 7, -1.5, -2.5, 2*1. 0, 8, -0.5, -2.5, 2*1. 0, 9, -3., -1., 2*2. 0, 10, -3.5, -1.5, 2*1. 0, 11, -2.5, -1.5, 2*1. 0, 12, -3.5, -0.5, 2*1. 0, 13, -2.5, -0.5, 2*1. 0, 14, 2*-1.5, 2*1. 0, 15, 2*-1., 2*2. 0, 16, -0.5, -1.5, 2*1. 0, 17, -1.5, -0.5, 2*1. 0, 18, 2*-0.5, 2*1. 0, 19, 1., -3., 2*2. 0, 20, 0.5, -3.5, 2*1. 0, 21, 2., -2., 2*4. 0, 22, 1.5, -3.5, 2*1. 0, 23, 0.5, -2.5, 2*1. 0, 24, 1.5, -2.5, 2*1. 0, 25, 3., -3., 2*2. 0, 26, 0.5, -1.5, 2*1. 0, 27, 1., -1., 2*2. 0, 28, 1.5, -1.5, 2*1. 0, 29, 0.5, -0.5, 2*1. 0, 30, 1.5, -0.5, 2*1. 0, 31, 2.5, -1.5, 2*1. 0, 32, 3., -1., 2*2. 0, 33, 3.5, -1.5, 2*1. 0, 34, 2.5, -0.5, 2*1. 0, 35, 3.5, -0.5, 2*1. 0, 36, -2., 2., 2*4. 0, 37, -3.5, 0.5, 2*1. 0, 38, -3., 1., 2*2. 0, 39, -2.5, 0.5, 2*1. 0, 40, -3.5, 1.5, 2*1. 0, 41, -2.5, 1.5, 2*1. 0, 42, -1.5, 0.5, 2*1. 0, 43, -1., 1., 2*2. 0, 44, -0.5, 0.5, 2*1. 0, 45, -1.5, 1.5, 2*1. 0, 46, -0.5, 1.5, 2*1. 0, 47, -3., 3., 2*2. 0, 48, -1.5, 2.5, 2*1. 0, 49, -1., 3., 2*2. 0, 50, -0.5, 2.5, 2*1. 0, 51, -1.5, 3.5, 2*1. 0, 52, -0.5, 3.5, 2*1. 0, 53, 2*0.5, 2*1. 0, 54, 2*1., 2*2. 0, 55, 2*2., 2*4. 0, 56, 1.5, 0.5, 2*1. 0, 57, 0.5, 1.5, 2*1. 0, 58, 2*1.5, 2*1. 0, 59, 3., 1., 2*2. 0, 60, 2.5, 0.5, 2*1. 0, 61, 3.5, 0.5, 2*1. 0, 62, 2.5, 1.5, 2*1. 0, 63, 3.5, 1.5, 2*1. 0, 64, 1., 3., 2*2. 0, 65, 0.5, 2.5, 2*1. 0, 66, 1.5, 2.5, 2*1. 0, 67, 0.5, 3.5, 2*1. 0, 68, 1.5, 3.5, 2*1. 0, 69, 2*3., 2*2.
Step 24. - Check solution.
do lb=1,lnblocks
if(coord(1,lb).eq..5.and.coord(2,lb).eq..5) then
do j=1,nyb+2*nguard
write(*,50) j,(unk(1,i,j,1,lb),i=1,nxb+2*nguard)
enddo
endif
enddo
gmake -f make_tutor your_tutorial
mpirun -np 1 your_tutorial/tutor
Explanation:
Your final lines of output should look like this if you are using
single precision,
1 2.6343 2.6343 2.6062 2.5514 2.4728 2.3744 2 2.6343 2.6343 2.6062 2.5514 2.4728 2.3744 3 2.6062 2.6062 2.5785 2.5247 2.4475 2.3508 4 2.5514 2.5514 2.5247 2.4727 2.3982 2.3049 5 2.4728 2.4728 2.4475 2.3982 2.3275 2.2391 6 2.3744 2.3744 2.3508 2.3049 2.2391 2.1567or like this with double precision (after modifying the format statement)
1 2.63432009 2.63432009 2.60617704 2.55138273 2.47279885 2.37442856 2 2.63432009 2.63432009 2.60617704 2.55138273 2.47279885 2.37442856 3 2.60617704 2.60617704 2.57852729 2.52469358 2.44748775 2.35084293 4 2.55138273 2.55138273 2.52469358 2.47273061 2.39820857 2.30492442 5 2.47279885 2.47279885 2.44748775 2.39820857 2.32753714 2.23907542 6 2.37442856 2.37442856 2.35084293 2.30492442 2.23907542 2.15665446The final solution is displayed here.
Step 25. - A Parallel Run
The code you have just developed will run in parallel.
Run the code using the command:
mpirun
-np X tutor where X
is the number of processors you wish to use.
Check that the last 6 lines of output are the same as before. If they are, then you have successfully completed this tutorial.
Step 26. - Adding Graphics