setup.c 4.01 KB
Newer Older
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
/* Setup routines for heat equation solver */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <assert.h>

#include "heat.h"
#include "pngwriter.h"

#define NSTEPS 500  // Default number of iteration steps

/* Initialize the heat equation solver */
void initialize(int argc, char *argv[], field *current,
                field *previous, int *nsteps)
{
    /*
     * Following combinations of command line arguments are possible:
     * No arguments:    use default field dimensions and number of time steps
     * One argument:    read initial field from a given file
     * Two arguments:   initial field from file and number of time steps
     * Three arguments: field dimensions (rows,cols) and number of time steps
     */


    int rows = 1024;             //!< Field dimensions with default values
    int cols = 1024;

    char input_file[64];        //!< Name of the optional input file

    int read_file = 0;

    *nsteps = NSTEPS;

    #pragma omp single
    {
        switch (argc) {
        case 1:
            /* Use default values */
            break;
        case 2:
            /* Read initial field from a file */
            strncpy(input_file, argv[1], 64);
            read_file = 1;
            break;
        case 3:
            /* Read initial field from a file */
            strncpy(input_file, argv[1], 64);
            read_file = 1;

            /* Number of time steps */
            *nsteps = atoi(argv[2]);
            break;
        case 4:
            /* Field dimensions */
            rows = atoi(argv[1]);
            cols = atoi(argv[2]);
            /* Number of time steps */
            *nsteps = atoi(argv[3]);
            break;
        default:
            printf("Unsupported number of command line arguments\n");
            exit(-1);
        }
    }

    if (read_file)
        #pragma omp single
        read_field(current, previous, input_file);
    else {
        #pragma omp single
        {
            set_field_dimensions(current, rows, cols);
            set_field_dimensions(previous, rows, cols);
        }
        generate_field(current);
        generate_field(previous);
    }
}

/* Generate initial temperature field.  Pattern is disc with a radius
 * of nx_full / 6 in the center of the grid.
 * Boundary conditions are (different) constant temperatures outside the grid */
void generate_field(field *temperature)
{
    int i, j;
    double radius;
    int dx, dy;

    /* Allocate the temperature array, note that
     * we have to allocate also the ghost layers */
    #pragma omp single
    {
        temperature->data =
            malloc_2d(temperature->nx + 2, temperature->ny + 2);
    }

    /* Radius of the source disc */
    radius = temperature->nx / 6.0;
    #pragma omp for private(i, j, dx, dy)
    for (i = 0; i < temperature->nx + 2; i++) {
        for (j = 0; j < temperature->ny + 2; j++) {
            /* Distances of point i, j from the origin */
            dx = i - temperature->nx / 2 + 1;
            dy = j - temperature->ny / 2 + 1;
            if (dx * dx + dy * dy < radius * radius) {
                temperature->data[i][j] = 5.0;
            } else {
                temperature->data[i][j] = 65.0;
            }
        }
    }

    /* Boundary conditions */
    #pragma omp for private(i)
    for (i = 0; i < temperature->nx + 2; i++) {
        temperature->data[i][0] = 20.0;
        temperature->data[i][temperature->ny + 1] = 70.0;
    }

    #pragma omp for private(j)
    for (j = 0; j < temperature->ny + 2; j++) {
        temperature->data[0][j] = 85.0;
        temperature->data[temperature->nx + 1][j] = 5.0;
    }

}

/* Set dimensions of the field. Note that the nx is the size of the first
 * dimension and ny the second. */
void set_field_dimensions(field *temperature, int nx, int ny)
{
    temperature->dx = DX;
    temperature->dy = DY;
    temperature->nx = nx;
    temperature->ny = ny;
}

/* Deallocate the 2D arrays of temperature fields */
void finalize(field *temperature1, field *temperature2)
{
    free_2d(temperature1->data);
    free_2d(temperature2->data);
}