#!/usr/bin/perl -w
use strict;
use Config;
use Cwd;
my $file = "heat_plate_mpi.pl";
local(*OUTF);
open(OUTF, ">$file") or die "Cannot open $file for writing: $!\n";
print OUTF $Config{startperl}, "\n\n";
print OUTF "use lib qw(", Cwd::cwd, "/../../blib/arch ",
Cwd::cwd, "/../../blib/lib);\n\n";
print "Writing $file\n";
while(<DATA>) { print OUTF $_ }
close(OUTF);
chmod(0755, $file);
__END__
$| = 1;
use Parallel::MPI qw(:all);
use vars qw($t $s $r $t1 $t2 $h $xl $yl $tol $h2 $aa $sum $global_sum
$nx $ny $i $j $ij $nproc $iter $last_iter $niter $n_left $npt $nxy
$nproc $my_rank $root @istatus @local_n @local_a);
MPI_Init();
$nproc = MPI_Comm_size(MPI_COMM_WORLD);
$my_rank = MPI_Comm_rank(MPI_COMM_WORLD);
$xl = 10.0;
$yl = 10.0;
$tol = 0.001;
# initial set up done by my_rank = 0 (i.e. master) process only
if($my_rank eq 0) {
print "there are $nproc processes\n";
print "enter number of iterations and h:";
($niter, $h) = split(' ', scalar(<STDIN>));
# determine the workload for each processor
$nx = $xl/$h + 1;
$ny = $nx;
print "nx=ny=$nx\n";
foreach(0..($nproc-1)) {
use integer;
$local_n[$_] = $nx/$nproc;
}
# if workload cannot be divided evenly
# to each process, assign to the
# first n_left processes with an
# additional piece of work
if(($nx % $nproc) != 0) {
use integer;
$n_left = $nx - $nproc*($nx/$nproc);
foreach(0..($n_left-1)) {
$local_n[$_]++;
}
$local_a[0] = 0.0;
foreach(1..($nproc-1)) {
$local_a[$_] = $local_a[$_-1] + $h*$local_n[i-1]
}
}
print "process #, local workload, start location\n";
foreach(0..($nproc-1)) {
print $_, "\t", $local_n[$_], "\t", $local_a[$_], "\n";
}
}
# now all the processes begin to work together
MPI_Bcast(\$h, 1,MPI_DOUBLE,0,MPI_COMM_WORLD);
MPI_Bcast(\$niter,1,MPI_INT, 0,MPI_COMM_WORLD);
MPI_Scatter(\@local_n,1,MPI_INT, \$npt,1,MPI_INT, 0,MPI_COMM_WORLD);
MPI_Scatter(\@local_a,1,MPI_DOUBLE,\$aa, 1,MPI_DOUBLE,0,MPI_COMM_WORLD);
# add one additional point to the end strips and
# add additional 2 points to the interior strips
# to facilitate overlap
$yl = 10.0;
$ny = $yl/$h + 1;
if($my_rank == 0 or $my_rank == ($nproc-1)) {
$nx = $npt + 1;
} else {
$nx = $npt + 2;
}
print "I am process $my_rank nx,ny= $nx,$ny\n";
print "[$my_rank] aa=$aa\n";
print "[$my_rank] h=$h\n";
print "[$my_rank] niter=$niter\n";
$h2 = $h*$h;
$nxy = $nx*$ny;
# source term
foreach $j (0..($ny-1)) {
foreach $i (0..($nx-1)) {
$t[$i][$j] = 20.0;
$s[$i][$j] = 0.0;
$r[$i][$j] = 0.0;
}
}
# fix the boundary conditions
#
# domain is divided into horizontal stripes
# along the west and east walls
#foreach $j (0..($ny-1)) {
# $t[0][$j] = 20.0;
# $t[$nx-1][$j] = 20.0;
#}
# along the south and north walls
# south boundary is owned by the first processor
if($my_rank == 0) {
foreach $i (0..($nx-1)) {
$t[$i][1] = 20.0;
}
}
# north boundary is owned by the last processor
if($my_rank == ($nproc-1)) {
foreach $i (0..($ny-1)) {
$xx = $i * $h;
if($xx >= 3.0 and $xx <= 7.0) {
$t[$nx-1][$i] = 100.0;
} else {
$t[$nx-1][$i] = 20.0;
}
}
}
MPI_Barrier(MPI_COMM_WORLD);
if($my_rank == 0) {
print "finish updating boundary conditions\n";
$t1 = MPI_Wtime()
}
foreach $iter (1..$niter) {
# send and receive interfacial values from nearest neighbours
if($my_rank != 0 and $my_rank != ($nproc-1)) {
# along the south side
@status = MPI_Sendrecv($t[1],$ny,MPI_DOUBLE,$my_rank-1,10,
$t[0],$ny,MPI_DOUBLE,$my_rank-1,20,
MPI_COMM_WORLD);
# along the north side
@status = MPI_Sendrecv($t[$nx-2],$ny,MPI_DOUBLE,$my_rank+1,20,
$t[$nx-1],$ny,MPI_DOUBLE,$my_rank+1,10,
MPI_COMM_WORLD);
}
if($my_rank == 0) {
# along the north side
$root = 1;
@status = MPI_Sendrecv($t[$nx-2],$ny,MPI_DOUBLE,$root,20,
$t[$nx-1],$ny,MPI_DOUBLE,$root,10,
MPI_COMM_WORLD);
}
if($my_rank == ($nproc-1)) {
# along the south side
$root = $nproc-2;
@status = MPI_Sendrecv($t[1],$ny,MPI_DOUBLE,$root,10,
$t[0],$ny,MPI_DOUBLE,$root,20,
MPI_COMM_WORLD);
}
$sum = 0.0;
foreach $j (1..($ny-2)) {
foreach $i (1..($nx-2)) {
$r[$i][$j] = $s[$i][$j] * $h2 - (
$t[$i+1][$j] + $t[$i][$j+1]
- 4.0 * $t[$i][$j] +
$t[$i-1][$j] + $t[$i][$j-1]);
$sum = $sum + $r[$i][$j] * $r[$i][$j];
}
}
# update temperatures
foreach $j (0..($ny-1)) {
foreach $i (0..($nx-1)) {
$t[$i][$j] = $t[$i][$j] - 0.25 * $r[$i][$j];
}
}
# obtain global sum
MPI_Allreduce(\$sum,\$global_sum,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
$global_sum = sqrt($global_sum);
if($my_rank == 0) {
if(($iter % 50) == 0) {
print "iter=$iter\nglobal_sum=$global_sum\n";
}
}
$last_iter = $iter;
last if $global_sum <= $tol;
}
if($my_rank == 0) {
$t2 = MPI_Wtime();
print "no. of iterations took = $last_iter\n";
print "final L2 norm of residual: $global_sum\n";
print "wall clock time in seconds= ",($t2-$t1), "\n";
}
#sleep($my_rank);
#if($my_rank == 0) {
# foreach $i (0..$nx-2) {
# print "[$my_rank]";
# foreach $j (0..$ny-1) {
# printf "%3d ", $t[$i][$j];
# }
# print "\n";
# }
#} elsif($my_rank == $nproc-1) {
# foreach $i (1..$nx-1) {
# print "[$my_rank]";
# foreach $j (0..$ny-1) {
# printf "%3d ", $t[$i][$j];
# }
# print "\n";
# }
#} else {
# foreach $i (1..$nx-2) {
# print "[$my_rank]";
# foreach $j (0..$ny-1) {
# printf "%3d ", $t[$i][$j];
# }
# print "\n";
# }
#}
#print "\n" if $my_rank == 1;
MPI_Finalize();