forked from moebiusband/NuSiF-Solver
		
	Add VTK output for sequential case
This commit is contained in:
		@@ -18,6 +18,56 @@
 | 
			
		||||
#include "timing.h"
 | 
			
		||||
#include "vtkWriter.h"
 | 
			
		||||
 | 
			
		||||
#ifndef _MPI
 | 
			
		||||
#define G(v, i, j, k) v[(k) * (imax + 2) * (jmax + 2) + (j) * (imax + 2) + (i)]
 | 
			
		||||
 | 
			
		||||
static void createBulkArrays(Solver* s, double* pg, double* ug, double* vg, double* wg)
 | 
			
		||||
{
 | 
			
		||||
    int imax = s->grid.imax;
 | 
			
		||||
    int jmax = s->grid.jmax;
 | 
			
		||||
    int kmax = s->grid.kmax;
 | 
			
		||||
    int idx  = 0;
 | 
			
		||||
 | 
			
		||||
    for (int k = 1; k < kmax + 1; k++) {
 | 
			
		||||
        for (int j = 1; j < jmax + 1; j++) {
 | 
			
		||||
            for (int i = 1; i < imax + 1; i++) {
 | 
			
		||||
                pg[idx++] = G(s->p, i, j, k);
 | 
			
		||||
            }
 | 
			
		||||
        }
 | 
			
		||||
    }
 | 
			
		||||
 | 
			
		||||
    idx = 0;
 | 
			
		||||
 | 
			
		||||
    for (int k = 1; k < kmax + 1; k++) {
 | 
			
		||||
        for (int j = 1; j < jmax + 1; j++) {
 | 
			
		||||
            for (int i = 1; i < imax + 1; i++) {
 | 
			
		||||
                ug[idx++] = (G(s->u, i, j, k) + G(s->u, i - 1, j, k)) / 2.0;
 | 
			
		||||
            }
 | 
			
		||||
        }
 | 
			
		||||
    }
 | 
			
		||||
 | 
			
		||||
    idx = 0;
 | 
			
		||||
 | 
			
		||||
    for (int k = 1; k < kmax + 1; k++) {
 | 
			
		||||
        for (int j = 1; j < jmax + 1; j++) {
 | 
			
		||||
            for (int i = 1; i < imax + 1; i++) {
 | 
			
		||||
                vg[idx++] = (G(s->v, i, j, k) + G(s->v, i, j - 1, k)) / 2.0;
 | 
			
		||||
            }
 | 
			
		||||
        }
 | 
			
		||||
    }
 | 
			
		||||
 | 
			
		||||
    idx = 0;
 | 
			
		||||
 | 
			
		||||
    for (int k = 1; k < kmax + 1; k++) {
 | 
			
		||||
        for (int j = 1; j < jmax + 1; j++) {
 | 
			
		||||
            for (int i = 1; i < imax + 1; i++) {
 | 
			
		||||
                wg[idx++] = (G(s->w, i, j, k) + G(s->w, i, j, k - 1)) / 2.0;
 | 
			
		||||
            }
 | 
			
		||||
        }
 | 
			
		||||
    }
 | 
			
		||||
}
 | 
			
		||||
#endif /* ifndef _MPI */
 | 
			
		||||
 | 
			
		||||
int main(int argc, char** argv)
 | 
			
		||||
{
 | 
			
		||||
    double timeStart, timeStop;
 | 
			
		||||
@@ -84,6 +134,7 @@ int main(int argc, char** argv)
 | 
			
		||||
        wg = allocate(64, bytesize);
 | 
			
		||||
    }
 | 
			
		||||
 | 
			
		||||
#ifdef _MPI
 | 
			
		||||
    commCollectResult(&s.comm,
 | 
			
		||||
        ug,
 | 
			
		||||
        vg,
 | 
			
		||||
@@ -96,12 +147,17 @@ int main(int argc, char** argv)
 | 
			
		||||
        s.grid.kmax,
 | 
			
		||||
        s.grid.jmax,
 | 
			
		||||
        s.grid.imax);
 | 
			
		||||
#else
 | 
			
		||||
    createBulkArrays(&s, pg, ug, vg, wg);
 | 
			
		||||
#endif /* ifdef _MPI */
 | 
			
		||||
 | 
			
		||||
    VtkOptions opts = { .grid = s.grid, .comm = s.comm };
 | 
			
		||||
    vtkOpen(&opts, s.problem);
 | 
			
		||||
    vtkScalar(&opts, "pressure", pg);
 | 
			
		||||
    vtkVector(&opts, "velocity", (VtkVector) { ug, vg, wg });
 | 
			
		||||
    vtkClose(&opts);
 | 
			
		||||
    if (commIsMaster(&s.comm)) {
 | 
			
		||||
        VtkOptions opts = { .grid = s.grid };
 | 
			
		||||
        vtkOpen(&opts, s.problem);
 | 
			
		||||
        vtkScalar(&opts, "pressure", pg);
 | 
			
		||||
        vtkVector(&opts, "velocity", (VtkVector) { ug, vg, wg });
 | 
			
		||||
        vtkClose(&opts);
 | 
			
		||||
    }
 | 
			
		||||
 | 
			
		||||
    commFinalize(&s.comm);
 | 
			
		||||
    return EXIT_SUCCESS;
 | 
			
		||||
 
 | 
			
		||||
@@ -55,16 +55,11 @@ void vtkOpen(VtkOptions* o, char* problem)
 | 
			
		||||
{
 | 
			
		||||
    char filename[50];
 | 
			
		||||
 | 
			
		||||
    if (o->mode == UNIX) {
 | 
			
		||||
        if (commIsMaster(&o->comm)) {
 | 
			
		||||
            snprintf(filename, 50, "%s-p%d.vtk", problem, o->comm.size);
 | 
			
		||||
            o->fh = fopen(filename, "w");
 | 
			
		||||
            writeHeader(o);
 | 
			
		||||
        }
 | 
			
		||||
    } else if (o->mode == MPI) {
 | 
			
		||||
    }
 | 
			
		||||
    snprintf(filename, 50, "%s.vtk", problem);
 | 
			
		||||
    o->fh = fopen(filename, "w");
 | 
			
		||||
    writeHeader(o);
 | 
			
		||||
 | 
			
		||||
    if (commIsMaster(&o->comm)) printf("Writing VTK output for %s\n", problem);
 | 
			
		||||
    printf("Writing VTK output for %s\n", problem);
 | 
			
		||||
}
 | 
			
		||||
 | 
			
		||||
static void writeScalar(VtkOptions* o, double* s)
 | 
			
		||||
@@ -101,17 +96,11 @@ static bool isInitialized(FILE* ptr)
 | 
			
		||||
 | 
			
		||||
void vtkScalar(VtkOptions* o, char* name, double* s)
 | 
			
		||||
{
 | 
			
		||||
    if (commIsMaster(&o->comm)) printf("Register scalar %s\n", name);
 | 
			
		||||
 | 
			
		||||
    if (o->mode == UNIX) {
 | 
			
		||||
        if (commIsMaster(&o->comm)) {
 | 
			
		||||
            if (!isInitialized(o->fh)) return;
 | 
			
		||||
            fprintf(o->fh, "SCALARS %s double 1\n", name);
 | 
			
		||||
            fprintf(o->fh, "LOOKUP_TABLE default\n");
 | 
			
		||||
            writeScalar(o, s);
 | 
			
		||||
        }
 | 
			
		||||
    } else if (o->mode == MPI) {
 | 
			
		||||
    }
 | 
			
		||||
    printf("Register scalar %s\n", name);
 | 
			
		||||
    if (!isInitialized(o->fh)) return;
 | 
			
		||||
    fprintf(o->fh, "SCALARS %s double 1\n", name);
 | 
			
		||||
    fprintf(o->fh, "LOOKUP_TABLE default\n");
 | 
			
		||||
    writeScalar(o, s);
 | 
			
		||||
}
 | 
			
		||||
 | 
			
		||||
static void writeVector(VtkOptions* o, VtkVector vec)
 | 
			
		||||
@@ -145,25 +134,14 @@ static void writeVector(VtkOptions* o, VtkVector vec)
 | 
			
		||||
 | 
			
		||||
void vtkVector(VtkOptions* o, char* name, VtkVector vec)
 | 
			
		||||
{
 | 
			
		||||
    if (commIsMaster(&o->comm)) printf("Register vector %s\n", name);
 | 
			
		||||
 | 
			
		||||
    if (o->mode == UNIX) {
 | 
			
		||||
        if (commIsMaster(&o->comm)) {
 | 
			
		||||
            if (!isInitialized(o->fh)) return;
 | 
			
		||||
            fprintf(o->fh, "VECTORS %s double\n", name);
 | 
			
		||||
            writeVector(o, vec);
 | 
			
		||||
        }
 | 
			
		||||
    } else if (o->mode == MPI) {
 | 
			
		||||
    }
 | 
			
		||||
    printf("Register vector %s\n", name);
 | 
			
		||||
    if (!isInitialized(o->fh)) return;
 | 
			
		||||
    fprintf(o->fh, "VECTORS %s double\n", name);
 | 
			
		||||
    writeVector(o, vec);
 | 
			
		||||
}
 | 
			
		||||
 | 
			
		||||
void vtkClose(VtkOptions* o)
 | 
			
		||||
{
 | 
			
		||||
    if (o->mode == UNIX) {
 | 
			
		||||
        if (commIsMaster(&o->comm)) {
 | 
			
		||||
            fclose(o->fh);
 | 
			
		||||
            o->fh = NULL;
 | 
			
		||||
        }
 | 
			
		||||
    } else if (o->mode == MPI) {
 | 
			
		||||
    }
 | 
			
		||||
    fclose(o->fh);
 | 
			
		||||
    o->fh = NULL;
 | 
			
		||||
}
 | 
			
		||||
 
 | 
			
		||||
@@ -8,18 +8,14 @@
 | 
			
		||||
#define __VTKWRITER_H_
 | 
			
		||||
#include <stdio.h>
 | 
			
		||||
 | 
			
		||||
#include "comm.h"
 | 
			
		||||
#include "grid.h"
 | 
			
		||||
 | 
			
		||||
typedef enum VtkFormat { ASCII = 0, BINARY } VtkFormat;
 | 
			
		||||
typedef enum VtkMode { UNIX = 0, MPI } VtkMode;
 | 
			
		||||
 | 
			
		||||
typedef struct VtkOptions {
 | 
			
		||||
    VtkFormat fmt;
 | 
			
		||||
    VtkMode mode;
 | 
			
		||||
    Grid grid;
 | 
			
		||||
    FILE* fh;
 | 
			
		||||
    Comm comm;
 | 
			
		||||
} VtkOptions;
 | 
			
		||||
 | 
			
		||||
typedef struct VtkVector {
 | 
			
		||||
 
 | 
			
		||||
		Reference in New Issue
	
	Block a user