ex7: init
This commit is contained in:
@@ -0,0 +1,4 @@
|
||||
data/
|
||||
sequential
|
||||
parallel
|
||||
*.png
|
||||
@@ -0,0 +1,30 @@
|
||||
CC=gcc
|
||||
PARALLEL_CC=nvcc
|
||||
CFLAGS+= -std=c99 -O2 -Wall -Wextra
|
||||
LDLIBS+= -lm
|
||||
SEQUENTIAL_SRC_FILES=wave_2d_sequential.c
|
||||
PARALLEL_SRC_FILES=wave_2d_parallel.cu
|
||||
IMAGES=$(shell find data -type f | sed s/\\.dat/.png/g | sed s/data/images/g )
|
||||
.PHONY: all clean dirs plot movie
|
||||
all: dirs ${TARGETS}
|
||||
dirs:
|
||||
mkdir -p data images
|
||||
sequential: ${SEQUENTIAL_SRC_FILES}
|
||||
$(CC) $^ $(CFLAGS) -o $@ $(LDLIBS)
|
||||
parallel: ${PARALLEL_SRC_FILES}
|
||||
$(PARALLEL_CC) $^ -O2 -fmad=false -o $@ $(LDLIBS)
|
||||
plot: ${IMAGES}
|
||||
images/%.png: data/%.dat
|
||||
./plot_image.sh $<
|
||||
movie: ${IMAGES}
|
||||
ffmpeg -y -an -i images/%5d.png -vcodec libx264 -pix_fmt yuv420p -profile:v baseline -level 3 -r 12 wave.mp4
|
||||
check: dirs sequential parallel
|
||||
mkdir -p data_sequential
|
||||
./sequential
|
||||
cp -rf ./data/* ./data_sequential
|
||||
./parallel
|
||||
python3 compare.py data_sequential/00000.dat data/00000.dat
|
||||
python3 compare.py data_sequential/00075.dat data/00075.dat
|
||||
rm -rf data_sequential
|
||||
clean:
|
||||
-rm -fr sequential parallel data images wave.mp4
|
||||
@@ -0,0 +1,7 @@
|
||||
* make : Creates the output folders
|
||||
* make parallel : Builds the parallel version
|
||||
* make sequential : builds the sequential version
|
||||
* ./parallel : Runs the parallel version and fills the 'data/' directory with stored time steps
|
||||
* make plot : converts saved time steps to png files under 'images/', using gnuplot. Runs faster if launched with e.g. 4 threads (make -j4 plot).
|
||||
* make movie : converts collection of png files under 'images' into an mp4 movie file, using ffmpeg
|
||||
* make check : builds both executeables and compares their output
|
||||
@@ -0,0 +1,48 @@
|
||||
import sys
|
||||
|
||||
import numpy as np
|
||||
|
||||
|
||||
class bcolors:
|
||||
OKGREEN = "\033[92m"
|
||||
FAIL = "\033[91m"
|
||||
ENDC = "\033[0m"
|
||||
|
||||
|
||||
def main() -> None:
|
||||
error_margin = 1e-4
|
||||
|
||||
if len(sys.argv) != 3:
|
||||
print("Usage: python plot_difference.py file1.dat file2.dat")
|
||||
sys.exit(1)
|
||||
|
||||
file1, file2 = sys.argv[1], sys.argv[2]
|
||||
|
||||
M = 128
|
||||
N = 128
|
||||
|
||||
data1 = load_data(file1, M, N)
|
||||
data2 = load_data(file2, M, N)
|
||||
|
||||
# Compute difference
|
||||
diff = data2 - data1
|
||||
|
||||
if (diff > error_margin).any():
|
||||
print(
|
||||
f"\n{bcolors.FAIL}Data files {file1} and {file2} differ by more than {error_margin}{bcolors.ENDC}\n"
|
||||
)
|
||||
else:
|
||||
print(
|
||||
f"\n{bcolors.OKGREEN}Data files {file1} and {file2} are identical within the margin of {error_margin}{bcolors.ENDC}\n"
|
||||
)
|
||||
|
||||
|
||||
def load_data(filename: str, M: int, N: int) -> np.ndarray:
|
||||
data = np.fromfile(filename, dtype=np.float64)
|
||||
if data.size != M * N:
|
||||
raise ValueError(f"File {filename} does not contain M*N={M*N} entries")
|
||||
return data.reshape((M, N))
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
@@ -0,0 +1,42 @@
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
import sys
|
||||
|
||||
|
||||
def main() -> None:
|
||||
if len(sys.argv) != 3:
|
||||
print("Usage: python plot_difference.py file1.dat file2.dat")
|
||||
sys.exit(1)
|
||||
|
||||
file1, file2 = sys.argv[1], sys.argv[2]
|
||||
|
||||
M = 128
|
||||
N = 128
|
||||
|
||||
data1 = load_data(file1, M, N)
|
||||
data2 = load_data(file2, M, N)
|
||||
|
||||
# Compute difference
|
||||
diff = data2 - data1
|
||||
|
||||
plt.figure(figsize=(6, 5))
|
||||
im = plt.imshow(diff, origin='lower', cmap='bwr', interpolation='nearest')
|
||||
plt.colorbar(im, label="Difference")
|
||||
plt.title(f"Difference between {file2} and {file1}")
|
||||
plt.xlabel("Column")
|
||||
plt.ylabel("Row")
|
||||
plt.tight_layout()
|
||||
|
||||
# Save figure to a file
|
||||
plt.savefig("difference_plot.png", dpi=150)
|
||||
print("Saved difference plot to difference_plot.png")
|
||||
|
||||
|
||||
def load_data(filename:str, M:int, N:int) -> np.ndarray:
|
||||
data = np.fromfile(filename, dtype=np.float64)
|
||||
if data.size != M * N:
|
||||
raise ValueError(f"File {filename} does not contain M*N={M*N} entries")
|
||||
return data.reshape((M, N))
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
Executable
+100
@@ -0,0 +1,100 @@
|
||||
#! /usr/bin/env bash
|
||||
|
||||
help()
|
||||
{
|
||||
echo
|
||||
echo "Plot 2D Wave Equation"
|
||||
echo
|
||||
echo "Syntax"
|
||||
echo "--------------------------------------------------------"
|
||||
echo "./plot_results.sh [-m|n|h] [data_folder] "
|
||||
echo
|
||||
echo "Option Description Arguments Default"
|
||||
echo "--------------------------------------------------------"
|
||||
echo "m y size Optional 128 "
|
||||
echo "n x size Optional 128 "
|
||||
echo "h Help None "
|
||||
echo
|
||||
echo "Example"
|
||||
echo "--------------------------------------------------------"
|
||||
echo "./plot_solution.sh -m 128 -n 128"
|
||||
echo
|
||||
}
|
||||
|
||||
#-----------------------------------------------------------------
|
||||
set -e
|
||||
|
||||
M=128
|
||||
N=128
|
||||
|
||||
# Check if the data folder is provided
|
||||
if [ $# -lt 1 ]; then
|
||||
echo "Error: No data folder provided."
|
||||
help
|
||||
exit 1
|
||||
fi
|
||||
|
||||
# Parse options and arguments
|
||||
while getopts ":m:n:h" opt; do
|
||||
case $opt in
|
||||
m)
|
||||
M=$OPTARG;;
|
||||
n)
|
||||
N=$OPTARG;;
|
||||
h)
|
||||
help
|
||||
exit;;
|
||||
\?)
|
||||
echo "Invalid option"
|
||||
help
|
||||
exit;;
|
||||
esac
|
||||
done
|
||||
|
||||
# Shift parsed options so that the remaining arguments start at $1
|
||||
shift $((OPTIND - 1))
|
||||
|
||||
# Ensure that the data folder is provided and exists
|
||||
DATAFOLDER=./data
|
||||
if [ ! -d "$DATAFOLDER" ]; then
|
||||
echo "Error: Data folder $DATAFOLDER does not exist."
|
||||
exit 1
|
||||
fi
|
||||
|
||||
#-----------------------------------------------------------------
|
||||
# Set up the size of the grid based on the options passed
|
||||
SIZE_M=`echo $M | bc`
|
||||
SIZE_N=`echo $N | bc`
|
||||
|
||||
# Ensure the output directory exists
|
||||
mkdir -p images
|
||||
|
||||
# Loop through all .dat files in the data folder
|
||||
for DATAFILE in "$DATAFOLDER"/*.dat; do
|
||||
# Skip if no .dat files are found
|
||||
if [ ! -f "$DATAFILE" ]; then
|
||||
echo "No .dat files found in the folder."
|
||||
exit 1
|
||||
fi
|
||||
|
||||
# Create the corresponding output image file name
|
||||
IMAGEFILE=`echo $DATAFILE | sed 's/dat$/png/' | sed 's/data/images/'`
|
||||
|
||||
# Run the gnuplot command to create the plot in the background
|
||||
(
|
||||
cat <<END_OF_SCRIPT | gnuplot -
|
||||
set term png
|
||||
set output "$IMAGEFILE"
|
||||
set zrange[-1:1]
|
||||
splot "$DATAFILE" binary array=${SIZE_M}x${SIZE_N} format='%double' with pm3d
|
||||
END_OF_SCRIPT
|
||||
|
||||
echo "Plot saved to $IMAGEFILE"
|
||||
) & # Run in the background
|
||||
|
||||
done
|
||||
|
||||
# Wait for all background processes to finish
|
||||
wait
|
||||
|
||||
echo "All plots have been generated."
|
||||
@@ -0,0 +1,258 @@
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
#include <stdint.h>
|
||||
#include <math.h>
|
||||
#include <sys/time.h>
|
||||
#include <sys/stat.h>
|
||||
#include <errno.h>
|
||||
#include <inttypes.h>
|
||||
|
||||
// TASK: T1
|
||||
// Include the cooperative groups library
|
||||
// BEGIN: T1
|
||||
;
|
||||
// END: T1
|
||||
|
||||
|
||||
// Convert 'struct timeval' into seconds in double prec. floating point
|
||||
#define WALLTIME(t) ((double)(t).tv_sec + 1e-6 * (double)(t).tv_usec)
|
||||
|
||||
|
||||
// Option to change numerical precision
|
||||
typedef int64_t int_t;
|
||||
typedef double real_t;
|
||||
|
||||
// TASK: T1b
|
||||
// Variables needed for implementation
|
||||
// BEGIN: T1b
|
||||
|
||||
// Simulation parameters: size, step count, and how often to save the state
|
||||
int_t
|
||||
N = 128,
|
||||
M = 128,
|
||||
max_iteration = 1000000,
|
||||
snapshot_freq = 1000;
|
||||
|
||||
// Wave equation parameters, time step is derived from the space step
|
||||
const real_t
|
||||
c = 1.0,
|
||||
dx = 1.0,
|
||||
dy = 1.0;
|
||||
real_t
|
||||
dt;
|
||||
|
||||
// Buffers for three time steps, indexed with 2 ghost points for the boundary
|
||||
real_t
|
||||
*buffers[3] = { NULL, NULL, NULL };
|
||||
|
||||
#define U_prv(i,j) buffers[0][((i)+1)*(N+2)+(j)+1]
|
||||
#define U(i,j) buffers[1][((i)+1)*(N+2)+(j)+1]
|
||||
#define U_nxt(i,j) buffers[2][((i)+1)*(N+2)+(j)+1]
|
||||
// END: T1b
|
||||
|
||||
#define cudaErrorCheck(ans) { gpuAssert((ans), __FILE__, __LINE__); }
|
||||
inline void gpuAssert(cudaError_t code, const char *file, int line, bool abort=true)
|
||||
{
|
||||
if (code != cudaSuccess) {
|
||||
fprintf(stderr,"GPUassert: %s %s %d\n", cudaGetErrorString(code), file, line);
|
||||
if (abort) exit(code);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// Rotate the time step buffers.
|
||||
void move_buffer_window ( void )
|
||||
{
|
||||
real_t *temp = buffers[0];
|
||||
buffers[0] = buffers[1];
|
||||
buffers[1] = buffers[2];
|
||||
buffers[2] = temp;
|
||||
}
|
||||
|
||||
|
||||
// Save the present time step in a numbered file under 'data/'
|
||||
void domain_save ( int_t step )
|
||||
{
|
||||
char filename[256];
|
||||
|
||||
// Ensure output directory exists (ignore error if it already exists)
|
||||
if (mkdir("data", 0755) != 0 && errno != EEXIST) {
|
||||
perror("mkdir data");
|
||||
exit(EXIT_FAILURE);
|
||||
}
|
||||
|
||||
snprintf(filename, sizeof(filename), "data/%05" PRId64 ".dat", step);
|
||||
|
||||
FILE *out = fopen(filename, "wb");
|
||||
if (out == NULL) {
|
||||
perror("fopen output file");
|
||||
fprintf(stderr, "Failed to open '%s' for writing.\n", filename);
|
||||
exit(EXIT_FAILURE);
|
||||
}
|
||||
|
||||
for ( int_t i = 0; i < M; ++i ) {
|
||||
size_t written = fwrite ( &U(i,0), sizeof(real_t), (size_t)N, out );
|
||||
if ( written != (size_t)N ) {
|
||||
perror("fwrite");
|
||||
fclose(out);
|
||||
exit(EXIT_FAILURE);
|
||||
}
|
||||
}
|
||||
|
||||
if ( fclose(out) != 0 ) {
|
||||
perror("fclose");
|
||||
exit(EXIT_FAILURE);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// TASK: T4
|
||||
// Get rid of all the memory allocations
|
||||
void domain_finalize ( void )
|
||||
{
|
||||
// BEGIN: T4
|
||||
free ( buffers[0] );
|
||||
free ( buffers[1] );
|
||||
free ( buffers[2] );
|
||||
// END: T4
|
||||
}
|
||||
|
||||
|
||||
// TASK: T6
|
||||
// Neumann (reflective) boundary condition
|
||||
// BEGIN: T6
|
||||
void boundary_condition ( void )
|
||||
{
|
||||
for ( int_t i=0; i<M; i++ )
|
||||
{
|
||||
U(i,-1) = U(i,1);
|
||||
U(i,N) = U(i,N-2);
|
||||
}
|
||||
for ( int_t j=0; j<N; j++ )
|
||||
{
|
||||
U(-1,j) = U(1,j);
|
||||
U(M,j) = U(M-2,j);
|
||||
}
|
||||
}
|
||||
// END: T6
|
||||
|
||||
|
||||
// TASK: T5
|
||||
// Integration formula
|
||||
// BEGIN: T5
|
||||
void time_step ( void )
|
||||
{
|
||||
for ( int_t i=0; i<M; i++ )
|
||||
{
|
||||
for ( int_t j=0; j<N; j++ )
|
||||
{
|
||||
U_nxt(i,j) = -U_prv(i,j) + 2.0*U(i,j)
|
||||
+ (dt*dt*c*c)/(dx*dy) * (
|
||||
U(i-1,j)+U(i+1,j)+U(i,j-1)+U(i,j+1)-4.0*U(i,j)
|
||||
);
|
||||
}
|
||||
}
|
||||
}
|
||||
// END: T5
|
||||
|
||||
|
||||
// TASK: T7
|
||||
// Main time integration.
|
||||
void simulate( void )
|
||||
{
|
||||
// BEGIN: T7
|
||||
// Go through each time step
|
||||
for ( int_t iteration=0; iteration<=max_iteration; iteration++ )
|
||||
{
|
||||
if ( (iteration % snapshot_freq)==0 )
|
||||
{
|
||||
domain_save ( iteration / snapshot_freq );
|
||||
}
|
||||
|
||||
// Derive step t+1 from steps t and t-1
|
||||
boundary_condition();
|
||||
time_step();
|
||||
|
||||
// Rotate the time step buffers
|
||||
move_buffer_window();
|
||||
}
|
||||
// END: T7
|
||||
}
|
||||
|
||||
|
||||
// TASK: T8
|
||||
// GPU occupancy
|
||||
void occupancy( void )
|
||||
{
|
||||
// BEGIN: T8
|
||||
;
|
||||
// END: T8
|
||||
}
|
||||
|
||||
|
||||
// TASK: T2
|
||||
// Make sure at least one CUDA-capable device exists
|
||||
static bool init_cuda()
|
||||
{
|
||||
// BEGIN: T2
|
||||
return true;
|
||||
// END: T2
|
||||
}
|
||||
|
||||
|
||||
// TASK: T3
|
||||
// Set up our three buffers, and fill two with an initial perturbation
|
||||
// Function to determine occupancy and optimal configuration
|
||||
|
||||
|
||||
void domain_initialize ( void )
|
||||
{
|
||||
// BEGIN: T3
|
||||
bool locate_cuda = init_cuda();
|
||||
if (!locate_cuda)
|
||||
{
|
||||
exit( EXIT_FAILURE );
|
||||
}
|
||||
|
||||
buffers[0] = (real_t *) malloc ( (M+2)*(N+2)*sizeof(real_t) );
|
||||
buffers[1] = (real_t *) malloc ( (M+2)*(N+2)*sizeof(real_t) );
|
||||
buffers[2] = (real_t *) malloc ( (M+2)*(N+2)*sizeof(real_t) );
|
||||
|
||||
for ( int_t i=0; i<M; i++ )
|
||||
{
|
||||
for ( int_t j=0; j<N; j++ )
|
||||
{
|
||||
// Calculate delta (radial distance) adjusted for M x N grid
|
||||
real_t delta = sqrt ( ((i - M/2.0) * (i - M/2.0)) / (real_t)M +
|
||||
((j - N/2.0) * (j - N/2.0)) / (real_t)N );
|
||||
U_prv(i,j) = U(i,j) = exp ( -4.0*delta*delta );
|
||||
}
|
||||
}
|
||||
|
||||
// Set the time step for 2D case
|
||||
dt = dx*dy / (c * sqrt (dx*dx+dy*dy));
|
||||
}
|
||||
|
||||
|
||||
int main ( void )
|
||||
{
|
||||
// Set up the initial state of the domain
|
||||
domain_initialize();
|
||||
|
||||
struct timeval t_start, t_end;
|
||||
|
||||
gettimeofday ( &t_start, NULL );
|
||||
simulate();
|
||||
gettimeofday ( &t_end, NULL );
|
||||
|
||||
printf ( "Total elapsed time: %lf seconds\n",
|
||||
WALLTIME(t_end) - WALLTIME(t_start)
|
||||
);
|
||||
|
||||
occupancy();
|
||||
|
||||
// Clean up and shut down
|
||||
domain_finalize();
|
||||
exit ( EXIT_SUCCESS );
|
||||
}
|
||||
@@ -0,0 +1,193 @@
|
||||
#define _XOPEN_SOURCE 600
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
#include <stdint.h>
|
||||
#include <math.h>
|
||||
#include <sys/time.h>
|
||||
#include <sys/stat.h>
|
||||
#include <errno.h>
|
||||
#include <inttypes.h>
|
||||
|
||||
|
||||
// Convert 'struct timeval' into seconds in double prec. floating point
|
||||
#define WALLTIME(t) ((double)(t).tv_sec + 1e-6 * (double)(t).tv_usec)
|
||||
|
||||
// Option to change numerical precision
|
||||
typedef int64_t int_t;
|
||||
typedef double real_t;
|
||||
|
||||
// Simulation parameters: size, step count, and how often to save the state
|
||||
int_t
|
||||
N = 128,
|
||||
M = 128,
|
||||
max_iteration = 1000000,
|
||||
snapshot_freq = 1000;
|
||||
|
||||
// Wave equation parameters, time step is derived from the space step
|
||||
const real_t
|
||||
c = 1.0,
|
||||
dx = 1.0,
|
||||
dy = 1.0;
|
||||
real_t
|
||||
dt;
|
||||
|
||||
// Buffers for three time steps, indexed with 2 ghost points for the boundary
|
||||
real_t
|
||||
*buffers[3] = { NULL, NULL, NULL };
|
||||
|
||||
#define U_prv(i,j) buffers[0][((i)+1)*(N+2)+(j)+1]
|
||||
#define U(i,j) buffers[1][((i)+1)*(N+2)+(j)+1]
|
||||
#define U_nxt(i,j) buffers[2][((i)+1)*(N+2)+(j)+1]
|
||||
|
||||
|
||||
// Rotate the time step buffers.
|
||||
void move_buffer_window ( void )
|
||||
{
|
||||
real_t *temp = buffers[0];
|
||||
buffers[0] = buffers[1];
|
||||
buffers[1] = buffers[2];
|
||||
buffers[2] = temp;
|
||||
}
|
||||
|
||||
|
||||
// Save the present time step in a numbered file under 'data/'
|
||||
void domain_save ( int_t step )
|
||||
{
|
||||
char filename[256];
|
||||
|
||||
// Ensure output directory exists (ignore error if it already exists)
|
||||
if (mkdir("data", 0755) != 0 && errno != EEXIST) {
|
||||
perror("mkdir data");
|
||||
exit(EXIT_FAILURE);
|
||||
}
|
||||
|
||||
snprintf(filename, sizeof(filename), "data/%05" PRId64 ".dat", step);
|
||||
|
||||
FILE *out = fopen(filename, "wb");
|
||||
if (out == NULL) {
|
||||
perror("fopen output file");
|
||||
fprintf(stderr, "Failed to open '%s' for writing.\n", filename);
|
||||
exit(EXIT_FAILURE);
|
||||
}
|
||||
|
||||
for ( int_t i = 0; i < M; ++i ) {
|
||||
size_t written = fwrite ( &U(i,0), sizeof(real_t), (size_t)N, out );
|
||||
if ( written != (size_t)N ) {
|
||||
perror("fwrite");
|
||||
fclose(out);
|
||||
exit(EXIT_FAILURE);
|
||||
}
|
||||
}
|
||||
|
||||
if ( fclose(out) != 0 ) {
|
||||
perror("fclose");
|
||||
exit(EXIT_FAILURE);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// Set up our three buffers, and fill two with an initial perturbation
|
||||
void domain_initialize ( void )
|
||||
{
|
||||
buffers[0] = malloc ( (M+2)*(N+2)*sizeof(real_t) );
|
||||
buffers[1] = malloc ( (M+2)*(N+2)*sizeof(real_t) );
|
||||
buffers[2] = malloc ( (M+2)*(N+2)*sizeof(real_t) );
|
||||
|
||||
for ( int_t i=0; i<M; i++ )
|
||||
{
|
||||
for ( int_t j=0; j<N; j++ )
|
||||
{
|
||||
// Calculate delta (radial distance) adjusted for M x N grid
|
||||
real_t delta = sqrt ( ((i - M/2.0) * (i - M/2.0)) / (real_t)M +
|
||||
((j - N/2.0) * (j - N/2.0)) / (real_t)N );
|
||||
U_prv(i,j) = U(i,j) = exp ( -4.0*delta*delta );
|
||||
}
|
||||
}
|
||||
|
||||
// Set the time step for 2D case
|
||||
dt = dx*dy / (c * sqrt (dx*dx+dy*dy));
|
||||
}
|
||||
|
||||
|
||||
// Get rid of all the memory allocations
|
||||
void domain_finalize ( void )
|
||||
{
|
||||
free ( buffers[0] );
|
||||
free ( buffers[1] );
|
||||
free ( buffers[2] );
|
||||
}
|
||||
|
||||
|
||||
// Integration formula (Eq. 9 from the pdf document)
|
||||
void time_step ( void )
|
||||
{
|
||||
for ( int_t i=0; i<M; i++ )
|
||||
{
|
||||
for ( int_t j=0; j<N; j++ )
|
||||
{
|
||||
U_nxt(i,j) = -U_prv(i,j) + 2.0*U(i,j)
|
||||
+ (dt*dt*c*c)/(dx*dy) * (
|
||||
U(i-1,j)+U(i+1,j)+U(i,j-1)+U(i,j+1)-4.0*U(i,j)
|
||||
);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// Neumann (reflective) boundary condition
|
||||
void boundary_condition ( void )
|
||||
{
|
||||
for ( int_t i=0; i<M; i++ )
|
||||
{
|
||||
U(i,-1) = U(i,1);
|
||||
U(i,N) = U(i,N-2);
|
||||
}
|
||||
for ( int_t j=0; j<N; j++ )
|
||||
{
|
||||
U(-1,j) = U(1,j);
|
||||
U(M,j) = U(M-2,j);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// Main time integration.
|
||||
void simulate( void )
|
||||
{
|
||||
// Go through each time step
|
||||
for ( int_t iteration=0; iteration<=max_iteration; iteration++ )
|
||||
{
|
||||
if ( (iteration % snapshot_freq)==0 )
|
||||
{
|
||||
domain_save ( iteration / snapshot_freq );
|
||||
}
|
||||
|
||||
// Derive step t+1 from steps t and t-1
|
||||
boundary_condition();
|
||||
time_step();
|
||||
|
||||
// Rotate the time step buffers
|
||||
move_buffer_window();
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
int main ( void )
|
||||
{
|
||||
// Set up the initial state of the domain
|
||||
domain_initialize();
|
||||
|
||||
struct timeval t_start, t_end;
|
||||
|
||||
gettimeofday ( &t_start, NULL );
|
||||
simulate();
|
||||
gettimeofday ( &t_end, NULL );
|
||||
|
||||
printf ( "Total elapsed time: %lf seconds\n",
|
||||
WALLTIME(t_end) - WALLTIME(t_start)
|
||||
);
|
||||
|
||||
// Clean up and shut down
|
||||
domain_finalize();
|
||||
exit ( EXIT_SUCCESS );
|
||||
}
|
||||
Reference in New Issue
Block a user