forked from moebiusband/NuSiF-Solver
		
	Cleanup. Port MultiGrid to 3D-seq
This commit is contained in:
		@@ -1,61 +0,0 @@
 | 
			
		||||
/*
 | 
			
		||||
 * Copyright (C)  NHR@FAU, University Erlangen-Nuremberg.
 | 
			
		||||
 * All rights reserved.
 | 
			
		||||
 * Use of this source code is governed by a MIT-style
 | 
			
		||||
 * license that can be found in the LICENSE file.
 | 
			
		||||
 */
 | 
			
		||||
#ifdef __linux__
 | 
			
		||||
#ifdef _OPENMP
 | 
			
		||||
#include <pthread.h>
 | 
			
		||||
#include <sched.h>
 | 
			
		||||
#include <stdio.h>
 | 
			
		||||
#include <stdlib.h>
 | 
			
		||||
#include <sys/syscall.h>
 | 
			
		||||
#include <sys/types.h>
 | 
			
		||||
#include <unistd.h>
 | 
			
		||||
 | 
			
		||||
#define MAX_NUM_THREADS 128
 | 
			
		||||
#define gettid()        syscall(SYS_gettid)
 | 
			
		||||
 | 
			
		||||
static int getProcessorID(cpu_set_t* cpu_set)
 | 
			
		||||
{
 | 
			
		||||
    int processorId;
 | 
			
		||||
 | 
			
		||||
    for (processorId = 0; processorId < MAX_NUM_THREADS; processorId++) {
 | 
			
		||||
        if (CPU_ISSET(processorId, cpu_set)) {
 | 
			
		||||
            break;
 | 
			
		||||
        }
 | 
			
		||||
    }
 | 
			
		||||
    return processorId;
 | 
			
		||||
}
 | 
			
		||||
 | 
			
		||||
int affinity_getProcessorId()
 | 
			
		||||
{
 | 
			
		||||
    cpu_set_t cpu_set;
 | 
			
		||||
    CPU_ZERO(&cpu_set);
 | 
			
		||||
    sched_getaffinity(gettid(), sizeof(cpu_set_t), &cpu_set);
 | 
			
		||||
 | 
			
		||||
    return getProcessorID(&cpu_set);
 | 
			
		||||
}
 | 
			
		||||
 | 
			
		||||
void affinity_pinThread(int processorId)
 | 
			
		||||
{
 | 
			
		||||
    cpu_set_t cpuset;
 | 
			
		||||
    pthread_t thread;
 | 
			
		||||
 | 
			
		||||
    thread = pthread_self();
 | 
			
		||||
    CPU_ZERO(&cpuset);
 | 
			
		||||
    CPU_SET(processorId, &cpuset);
 | 
			
		||||
    pthread_setaffinity_np(thread, sizeof(cpu_set_t), &cpuset);
 | 
			
		||||
}
 | 
			
		||||
 | 
			
		||||
void affinity_pinProcess(int processorId)
 | 
			
		||||
{
 | 
			
		||||
    cpu_set_t cpuset;
 | 
			
		||||
 | 
			
		||||
    CPU_ZERO(&cpuset);
 | 
			
		||||
    CPU_SET(processorId, &cpuset);
 | 
			
		||||
    sched_setaffinity(0, sizeof(cpu_set_t), &cpuset);
 | 
			
		||||
}
 | 
			
		||||
#endif /*_OPENMP*/
 | 
			
		||||
#endif /*__linux__*/
 | 
			
		||||
@@ -1,14 +0,0 @@
 | 
			
		||||
/*
 | 
			
		||||
 * Copyright (C)  NHR@FAU, University Erlangen-Nuremberg.
 | 
			
		||||
 * All rights reserved.
 | 
			
		||||
 * Use of this source code is governed by a MIT-style
 | 
			
		||||
 * license that can be found in the LICENSE file.
 | 
			
		||||
 */
 | 
			
		||||
#ifndef AFFINITY_H
 | 
			
		||||
#define AFFINITY_H
 | 
			
		||||
 | 
			
		||||
extern int affinity_getProcessorId();
 | 
			
		||||
extern void affinity_pinProcess(int);
 | 
			
		||||
extern void affinity_pinThread(int);
 | 
			
		||||
 | 
			
		||||
#endif /*AFFINITY_H*/
 | 
			
		||||
@@ -24,7 +24,6 @@ void initParameter(Parameter* param)
 | 
			
		||||
    param->re      = 100.0;
 | 
			
		||||
    param->gamma   = 0.9;
 | 
			
		||||
    param->tau     = 0.5;
 | 
			
		||||
    param->rho     = 0.99;
 | 
			
		||||
    param->levels  = 5;
 | 
			
		||||
}
 | 
			
		||||
 | 
			
		||||
@@ -81,7 +80,6 @@ void readParameter(Parameter* param, const char* filename)
 | 
			
		||||
            PARSE_REAL(u_init);
 | 
			
		||||
            PARSE_REAL(v_init);
 | 
			
		||||
            PARSE_REAL(p_init);
 | 
			
		||||
            PARSE_REAL(rho);
 | 
			
		||||
        }
 | 
			
		||||
    }
 | 
			
		||||
 | 
			
		||||
@@ -112,6 +110,5 @@ void printParameter(Parameter* param)
 | 
			
		||||
    printf("\tepsilon (stopping tolerance) : %f\n", param->eps);
 | 
			
		||||
    printf("\tgamma (stopping tolerance) : %f\n", param->gamma);
 | 
			
		||||
    printf("\tomega (SOR relaxation): %f\n", param->omg);
 | 
			
		||||
    printf("\trho (SOR relaxation): %f\n", param->rho);
 | 
			
		||||
    printf("\tMultiGrid levels : %d\n", param->levels);
 | 
			
		||||
}
 | 
			
		||||
 
 | 
			
		||||
@@ -9,6 +9,6 @@
 | 
			
		||||
 | 
			
		||||
extern void initProgress(double);
 | 
			
		||||
extern void printProgress(double);
 | 
			
		||||
extern void stopProgress();
 | 
			
		||||
extern void stopProgress(void);
 | 
			
		||||
 | 
			
		||||
#endif
 | 
			
		||||
 
 | 
			
		||||
@@ -8,16 +8,15 @@
 | 
			
		||||
#include <stdlib.h>
 | 
			
		||||
 | 
			
		||||
#include "allocate.h"
 | 
			
		||||
#include "discretization.h"
 | 
			
		||||
#include "parameter.h"
 | 
			
		||||
#include "solver.h"
 | 
			
		||||
#include "util.h"
 | 
			
		||||
 | 
			
		||||
#define FINEST_LEVEL 0
 | 
			
		||||
#define S(i, j)      s[(j) * (imax + 2) + (i)]
 | 
			
		||||
#define E(i, j)      e[(j) * (imax + 2) + (i)]
 | 
			
		||||
#define R(i, j)      r[(j) * (imax + 2) + (i)]
 | 
			
		||||
#define OLD(i, j)    old[(j) * (imax + 2) + (i)]
 | 
			
		||||
#define FINEST_LEVEL   0
 | 
			
		||||
#define COARSEST_LEVEL (solver->levels - 1)
 | 
			
		||||
#define S(i, j)        s[(j) * (imax + 2) + (i)]
 | 
			
		||||
#define E(i, j)        e[(j) * (imax + 2) + (i)]
 | 
			
		||||
#define R(i, j)        r[(j) * (imax + 2) + (i)]
 | 
			
		||||
#define OLD(i, j)      old[(j) * (imax + 2) + (i)]
 | 
			
		||||
 | 
			
		||||
static void restrictMG(Solver* s, int level, int imax, int jmax)
 | 
			
		||||
{
 | 
			
		||||
@@ -51,6 +50,7 @@ static void prolongate(Solver* s, int level, int imax, int jmax)
 | 
			
		||||
static void correct(Solver* s, double* p, int level, int imax, int jmax)
 | 
			
		||||
{
 | 
			
		||||
    double* e = s->e[level];
 | 
			
		||||
 | 
			
		||||
    for (int j = 1; j < jmax + 1; ++j) {
 | 
			
		||||
        for (int i = 1; i < imax + 1; ++i) {
 | 
			
		||||
            P(i, j) += E(i, j);
 | 
			
		||||
@@ -106,37 +106,8 @@ static double smooth(Solver* s, double* p, double* rhs, int level, int imax, int
 | 
			
		||||
    return res;
 | 
			
		||||
}
 | 
			
		||||
 | 
			
		||||
void initSolver(Solver* s, Discretization* d, Parameter* p)
 | 
			
		||||
{
 | 
			
		||||
    s->eps     = p->eps;
 | 
			
		||||
    s->omega   = p->omg;
 | 
			
		||||
    s->itermax = p->itermax;
 | 
			
		||||
    s->rho     = p->rho;
 | 
			
		||||
    s->levels  = p->levels;
 | 
			
		||||
    s->grid    = &d->grid;
 | 
			
		||||
 | 
			
		||||
    int imax   = s->grid->imax;
 | 
			
		||||
    int jmax   = s->grid->jmax;
 | 
			
		||||
    int levels = s->levels;
 | 
			
		||||
    printf("Using Multigrid solver with %d levels\n", levels);
 | 
			
		||||
 | 
			
		||||
    s->r = malloc(levels * sizeof(double*));
 | 
			
		||||
    s->e = malloc(levels * sizeof(double*));
 | 
			
		||||
 | 
			
		||||
    size_t size = (imax + 2) * (jmax + 2) * sizeof(double);
 | 
			
		||||
 | 
			
		||||
    for (int j = 0; j < levels; j++) {
 | 
			
		||||
        s->r[j] = allocate(64, size);
 | 
			
		||||
        s->e[j] = allocate(64, size);
 | 
			
		||||
 | 
			
		||||
        for (int i = 0; i < (imax + 2) * (jmax + 2); i++) {
 | 
			
		||||
            s->r[j][i] = 0.0;
 | 
			
		||||
            s->e[j][i] = 0.0;
 | 
			
		||||
        }
 | 
			
		||||
    }
 | 
			
		||||
}
 | 
			
		||||
 | 
			
		||||
double multiGrid(Solver* solver, double* p, double* rhs, int level, int imax, int jmax)
 | 
			
		||||
static double multiGrid(
 | 
			
		||||
    Solver* solver, double* p, double* rhs, int level, int imax, int jmax)
 | 
			
		||||
{
 | 
			
		||||
    double res = 0.0;
 | 
			
		||||
 | 
			
		||||
@@ -161,7 +132,7 @@ double multiGrid(Solver* solver, double* p, double* rhs, int level, int imax, in
 | 
			
		||||
    // TODO: What if there is a rest?
 | 
			
		||||
    multiGrid(solver,
 | 
			
		||||
        solver->e[level + 1],
 | 
			
		||||
        solver->r[level],
 | 
			
		||||
        solver->r[level + 1],
 | 
			
		||||
        level + 1,
 | 
			
		||||
        imax / 2,
 | 
			
		||||
        jmax / 2);
 | 
			
		||||
@@ -182,6 +153,35 @@ double multiGrid(Solver* solver, double* p, double* rhs, int level, int imax, in
 | 
			
		||||
    return res;
 | 
			
		||||
}
 | 
			
		||||
 | 
			
		||||
void initSolver(Solver* s, Discretization* d, Parameter* p)
 | 
			
		||||
{
 | 
			
		||||
    s->eps     = p->eps;
 | 
			
		||||
    s->omega   = p->omg;
 | 
			
		||||
    s->itermax = p->itermax;
 | 
			
		||||
    s->levels  = p->levels;
 | 
			
		||||
    s->grid    = &d->grid;
 | 
			
		||||
 | 
			
		||||
    int imax   = s->grid->imax;
 | 
			
		||||
    int jmax   = s->grid->jmax;
 | 
			
		||||
    int levels = s->levels;
 | 
			
		||||
    printf("Using Multigrid solver with %d levels\n", levels);
 | 
			
		||||
 | 
			
		||||
    s->r = malloc(levels * sizeof(double*));
 | 
			
		||||
    s->e = malloc(levels * sizeof(double*));
 | 
			
		||||
 | 
			
		||||
    size_t size = (imax + 2) * (jmax + 2) * sizeof(double);
 | 
			
		||||
 | 
			
		||||
    for (int j = 0; j < levels; j++) {
 | 
			
		||||
        s->r[j] = allocate(64, size);
 | 
			
		||||
        s->e[j] = allocate(64, size);
 | 
			
		||||
 | 
			
		||||
        for (int i = 0; i < (imax + 2) * (jmax + 2); i++) {
 | 
			
		||||
            s->r[j][i] = 0.0;
 | 
			
		||||
            s->e[j][i] = 0.0;
 | 
			
		||||
        }
 | 
			
		||||
    }
 | 
			
		||||
}
 | 
			
		||||
 | 
			
		||||
void solve(Solver* s, double* p, double* rhs)
 | 
			
		||||
{
 | 
			
		||||
    double res = multiGrid(s, p, rhs, 0, s->grid->imax, s->grid->jmax);
 | 
			
		||||
 
 | 
			
		||||
@@ -4,7 +4,6 @@
 | 
			
		||||
 * Use of this source code is governed by a MIT style
 | 
			
		||||
 * license that can be found in the LICENSE file.
 | 
			
		||||
 */
 | 
			
		||||
#include "discretization.h"
 | 
			
		||||
#include "solver.h"
 | 
			
		||||
#include "util.h"
 | 
			
		||||
 | 
			
		||||
 
 | 
			
		||||
		Reference in New Issue
	
	Block a user