/************************************************************************ | |
* * | |
* Commonwealth Scientific and Industrial Research Organisation (CSIRO) * | |
* - Division of Information Technology (DIT) * | |
* - Division of Atmospheric Research (DAR) * | |
* * | |
* Shallow water weather model - Distributed Memory Version * | |
* * | |
* Finite difference model of shallow water equations based on :- * | |
* "The dynamics of finite difference models of the shallow water * | |
* equations" by R. Sadourney, JAS, 32, 1975. * | |
* Code from:- * | |
* "An introduction to three-dimensional climate modelling" * | |
* by Washington and Parkinson * | |
* * | |
* Programmers = David Abramson (DIT) rcoda@koel.co.rmit.oz * | |
* = Paul Whiting (DIT) rcopw@koel.co.rmit.oz * | |
* = Martin Dix (DAR) mrd@koel.co.rmit.oz * | |
* Language = BSD c using Argonne NL macros * | |
* O/S = Unix System V * | |
* H/W = Encore Multimax 320 * | |
* * | |
************************************************************************/ | |
#include <math.h> | |
#include <mpi.h> | |
#include <stdio.h> | |
#include "decs.h" | |
MPI_Datatype * setup_res(); | |
extern void tstep( | |
int m_, | |
int n_, | |
float alpha_, | |
int jstart, | |
int jend, | |
float pold[n][m], | |
float uold[n][m], | |
float vold[n][m], | |
float p[n][m], | |
float u[n][m], | |
float v[n][m], | |
float pnew[n][m], | |
float unew[n][m], | |
float vnew[n][m], | |
float dpdt[n][m], | |
float dudt[n][m], | |
float dvdt[n][m], | |
int firststep, | |
float tdt); | |
void | |
worker() | |
{ | |
int firststep, ncycle; | |
float tdt, time; | |
int i,j,ip,jp,jstart,jend; | |
int prv; | |
int nxt; | |
int msg_type; | |
int master_id, my_id; | |
int nprocs; | |
int nbytes; | |
float p[n][m]; /* Pressure (or free surface height) */ | |
float u[n][m]; /* Zonal wind */ | |
float v[n][m]; /* Meridional wind */ | |
float psi[n][m]; /* Velocity streamfunction */ | |
float pold[n][m]; | |
float uold[n][m]; | |
float vold[n][m]; | |
float pnew[n][m]; | |
float unew[n][m]; | |
float vnew[n][m]; | |
float dpdt[n][m]; | |
float dudt[n][m]; /* Time tendency of u */ | |
float dvdt[n][m]; | |
float cu[n][m]; /* Mass weighted u */ | |
float cv[n][m]; /* Mass weighted v */ | |
float h[n][m]; | |
float z[n][m]; /* Potential enstrophy */ | |
float dummy1[m]; | |
float dummy2[n]; | |
float fsdx = 4./dx; | |
float fsdy = 4./dy; | |
int worker[4]; | |
float p_start[m]; | |
float u_start[m]; | |
float v_start[m]; | |
float psi_start[m]; | |
float pold_start[m]; | |
float uold_start[m]; | |
float vold_start[m]; | |
MPI_Datatype * res_type; | |
MPI_Status status; | |
/* | |
initialise control variables | |
*/ | |
firststep = 1; | |
ncycle = 0; | |
tdt = dt; | |
time = 0.; | |
/* | |
set up environment for worker | |
*/ | |
MPI_Comm_size(MPI_COMM_WORLD, &nprocs); | |
MPI_Comm_rank(MPI_COMM_WORLD, &my_id); | |
MPI_Recv(&worker, 4, MPI_INT, MPI_ANY_SOURCE, START_SIGNAL, | |
MPI_COMM_WORLD, &status); | |
prv = worker[PREV]; | |
nxt = worker[NEXT]; | |
jstart = worker[JSTART]; | |
jend = worker[JEND]; | |
/* | |
receive initialisation packets from master | |
*/ | |
for (i = 0; i < n; i++){ | |
MPI_Recv(&p_start, m, MPI_FLOAT, 0, P_ROW, | |
MPI_COMM_WORLD, &status); | |
acopy_one_to_two(p_start, p, i); | |
MPI_Recv(&u_start, m, MPI_FLOAT, 0, U_ROW, | |
MPI_COMM_WORLD, &status); | |
acopy_one_to_two(u_start, u, i); | |
MPI_Recv(&v_start, m, MPI_FLOAT, 0, V_ROW, | |
MPI_COMM_WORLD, &status); | |
acopy_one_to_two(v_start, v, i); | |
MPI_Recv(&psi_start, m, MPI_FLOAT, 0, PSI_ROW, | |
MPI_COMM_WORLD, &status); | |
acopy_one_to_two(psi_start, psi, i); | |
MPI_Recv(&pold_start, m, MPI_FLOAT, 0, POLD_ROW, | |
MPI_COMM_WORLD, &status); | |
acopy_one_to_two(pold_start, pold, i); | |
MPI_Recv(&uold_start, m, MPI_FLOAT, 0, UOLD_ROW, | |
MPI_COMM_WORLD, &status); | |
acopy_one_to_two(uold_start, uold, i); | |
MPI_Recv(&vold_start, m, MPI_FLOAT, 0, VOLD_ROW, | |
MPI_COMM_WORLD, &status); | |
acopy_one_to_two(vold_start, vold, i); | |
} | |
while (ncycle < itmax) { | |
/* | |
loop over latitudes calculating U, V, z and h | |
do the block of latitudes from jstart to jend inclusive | |
*/ | |
calc_load(prv, nxt, my_id, jstart, jend, p, u,v); | |
calcuvzh(jstart, jend, p, u, v, cu, cv, h, z, fsdx, fsdy); | |
calc_unload(prv, nxt, my_id, jstart, jend, cv, z); | |
/* | |
Calculate time tendencies of p, u and v | |
*/ | |
time_load(prv, nxt, my_id, jstart, jend, cu, cv, h, z); | |
timetend(jstart, jend, dpdt, dudt, dvdt, cu, cv, h, z); | |
time_unload(prv, nxt, my_id, jstart, jend, dvdt); | |
if ((my_id == 1) && (ncycle%mprint==0)) { | |
diag(ncycle, time, p, u, v, h, z); | |
} | |
time += dt; | |
tstep(m, n, alpha, | |
jstart, jend, pold, uold, vold, p, u, v, pnew, | |
unew, vnew, dpdt, dudt, dvdt, firststep, tdt); | |
if ( firststep ) { | |
/* Double tdt because all future steps are leapfrog */ | |
firststep = 0; | |
tdt = tdt+tdt; | |
} | |
ncycle++; | |
} /* End of time step loop */ | |
/* | |
send local data structures (results) back to master | |
*/ | |
res_type = setup_res(); | |
send_updated_ds(res_type, jstart, jend, p, P_ROW, 0); | |
send_updated_ds(res_type, jstart, jend, u, U_ROW, 0); | |
send_updated_ds(res_type, jstart, jend, v, V_ROW, 0); | |
send_updated_ds(res_type, jstart, jend, h, H_ROW, 0); | |
send_updated_ds(res_type, jstart, jend, z, Z_ROW, 0); | |
MPI_Recv(&worker, 4, MPI_INT, 0, END_SIGNAL, | |
MPI_COMM_WORLD, &status); | |
if (debug & debug_call) { | |
printf("worker %d sent TIDY_UP to master\n", my_id); | |
printf("worker %d got END_SIGNAL from master\n", my_id); | |
} | |
} | |
send_updated_ds(res_type, jstart, jend, ds, indx, master_id) | |
MPI_Datatype * res_type; | |
int jstart; | |
int jend; | |
float ds[n][m]; | |
int indx; | |
int master_id; | |
{ | |
int j; | |
struct res res; | |
MPI_Request rq[2]; | |
MPI_Status stat[2]; | |
for (j = jstart; j <= jend; j++) { | |
acopy_two_to_one(ds, res.row, j); | |
res.indx = j; | |
MPI_Send(&res, 1, *res_type, master_id, indx, | |
MPI_COMM_WORLD); | |
} | |
} | |
/* | |
this procedure does all the message passing before the call to _calcuvzh_ | |
*/ | |
calc_load(prv,nxt,my_id,jstart,jend,p,u,v) | |
int prv; | |
int nxt; | |
int my_id; | |
int jstart; | |
int jend; | |
float p[n][m]; | |
float u[n][m]; | |
float v[n][m]; | |
{ | |
neighbour_send(prv, my_id, CALC1a, p, jstart); | |
neighbour_send(prv, my_id, CALC1b, u, jstart); | |
neighbour_send(prv, my_id, CALC1c, v, jstart); | |
neighbour_receive(nxt, my_id, CALC1a, p, (jend+1) % n); | |
neighbour_receive(nxt, my_id, CALC1b, u, (jend+1) % n); | |
neighbour_receive(nxt, my_id, CALC1c, v, (jend+1) % n); | |
} | |
/* | |
this procedure does all the message passing after the call to _calcuvzh_ | |
*/ | |
calc_unload(prv,nxt,my_id,jstart,jend,cv,z) | |
int prv; | |
int nxt; | |
int my_id; | |
int jstart; | |
int jend; | |
float cv[n][m]; | |
float z[n][m]; | |
{ | |
neighbour_send(nxt, my_id, CALC2a, cv, (jend+1) % n); | |
neighbour_send(nxt, my_id, CALC2b, z, (jend+1) % n); | |
neighbour_receive(prv, my_id, CALC2a, cv, jstart); | |
neighbour_receive(prv, my_id, CALC2b, z, jstart); | |
} | |
/* | |
this procedure does all the message passing before the call to _timetend_ | |
*/ | |
time_load(prv,nxt,my_id,jstart,jend,cu,cv,h,z) | |
int prv; | |
int nxt; | |
int my_id; | |
int jstart; | |
int jend; | |
float cu[n][m]; | |
float cv[n][m]; | |
float h[n][m]; | |
float z[n][m]; | |
{ | |
neighbour_send(prv, my_id, TIME1a, cu, jstart); | |
neighbour_send(prv, my_id, TIME1b, cv, jstart); | |
neighbour_send(prv, my_id, TIME1c, h, jstart); | |
neighbour_send(prv, my_id, TIME1d, z, jstart); | |
neighbour_receive(nxt, my_id, TIME1a, cu, (jend+1) % n); | |
neighbour_receive(nxt, my_id, TIME1b, cv, (jend+1) % n); | |
neighbour_receive(nxt, my_id, TIME1c, h, (jend+1) % n); | |
neighbour_receive(nxt, my_id, TIME1d, z, (jend+1) % n); | |
} | |
/* | |
this procedure does all the message passing after the call to _timetend_ | |
*/ | |
time_unload(prv,nxt,tu_my_id,jstart,jend,dvdt) | |
int prv; | |
int nxt; | |
int tu_my_id; | |
int jstart; | |
int jend; | |
float dvdt[n][m]; | |
{ | |
neighbour_send(nxt, tu_my_id, TIME2, dvdt, (jend+1) % n); | |
neighbour_receive(prv, tu_my_id, TIME2, dvdt, jstart); | |
} | |
/* | |
this is a general purpose function for sending packets b/w workers | |
*/ | |
neighbour_send(ns_neighbour,ns_my_id,ns_rec_id,ns_ds,ns_edge) | |
int ns_neighbour; | |
int ns_my_id; | |
int ns_rec_id; | |
float ns_ds[n][m]; | |
int ns_edge; | |
{ | |
float ns_rec[m]; | |
MPI_Request rq; | |
acopy_two_to_one(ns_ds, ns_rec, ns_edge); | |
MPI_Isend(&ns_rec, m, MPI_FLOAT, ns_neighbour, ns_rec_id, | |
MPI_COMM_WORLD, &rq); | |
if (debug & debug_worker) | |
printf("worker %d sent packet %d to worker %d\n", ns_my_id, | |
ns_rec_id, ns_neighbour); | |
} | |
/* | |
this is a general purpose function for receiving packets b/w workers | |
*/ | |
neighbour_receive(nr_neighbour,nr_my_id,nr_rec_id,nr_ds,nr_edge) | |
int nr_neighbour; | |
int nr_my_id; | |
int nr_rec_id; | |
float nr_ds[n][m]; | |
int nr_edge; | |
{ | |
float nr_rec[m]; | |
MPI_Status status; | |
MPI_Recv(&nr_rec, m, MPI_FLOAT, nr_neighbour, nr_rec_id, | |
MPI_COMM_WORLD, &status); | |
acopy_one_to_two(nr_rec, nr_ds, nr_edge); | |
} |