Review particle tracer. Refactor.
This commit is contained in:
parent
ed929ab752
commit
6f3d9e73ef
@ -18,9 +18,10 @@ include $(MAKE_DIR)/include_$(TAG).mk
|
||||
INCLUDES += -I$(SRC_DIR) -I$(BUILD_DIR)
|
||||
|
||||
VPATH = $(SRC_DIR)
|
||||
SRC = $(wildcard $(SRC_DIR)/*.c)
|
||||
SRC = $(filter-out $(wildcard $(SRC_DIR)/*-*.c),$(wildcard $(SRC_DIR)/*.c))
|
||||
ASM = $(patsubst $(SRC_DIR)/%.c, $(BUILD_DIR)/%.s, $(SRC))
|
||||
OBJ = $(patsubst $(SRC_DIR)/%.c, $(BUILD_DIR)/%.o, $(SRC))
|
||||
OBJ += $(BUILD_DIR)/solver-$(SOLVER).o
|
||||
SOURCES = $(SRC) $(wildcard $(SRC_DIR)/*.h)
|
||||
CPPFLAGS := $(CPPFLAGS) $(DEFINES) $(OPTIONS) $(INCLUDES)
|
||||
|
||||
@ -49,6 +50,8 @@ distclean: clean
|
||||
@rm -f $(TARGET)
|
||||
@rm -f *.dat
|
||||
@rm -f *.png
|
||||
@rm -f ./vis_files/*.dat
|
||||
@rm -f ./vis_files/*.gif
|
||||
|
||||
info:
|
||||
$(info $(CFLAGS))
|
||||
|
@ -26,13 +26,13 @@ p_init 1.0 # initial value for pressure
|
||||
|
||||
xlength 30.0 # domain size in x-direction
|
||||
ylength 4.0 # domain size in y-direction
|
||||
imax 200 # number of interior cells in x-direction
|
||||
jmax 50 # number of interior cells in y-direction
|
||||
imax 256 # number of interior cells in x-direction
|
||||
jmax 64 # number of interior cells in y-direction
|
||||
|
||||
# Time Data:
|
||||
# ---------
|
||||
|
||||
te 60.0 # final time
|
||||
te 80.0 # final time
|
||||
dt 0.02 # time stepsize
|
||||
tau 0.5 # safety factor for time stepsize control (<0 constant delt)
|
||||
|
||||
@ -50,9 +50,9 @@ levels 5 # Multigrid levels
|
||||
# -----------------------
|
||||
|
||||
numberOfParticles 60
|
||||
startTime 100
|
||||
injectTimePeriod 2.0
|
||||
writeTimePeriod 0.5
|
||||
startTime 5.0
|
||||
injectTimePeriod 4.0
|
||||
writeTimePeriod 1.0
|
||||
|
||||
x1 1.0
|
||||
y1 0.0
|
||||
|
@ -1,12 +1,12 @@
|
||||
# Supported: GCC, CLANG, ICC
|
||||
TAG ?= CLANG
|
||||
ENABLE_OPENMP ?= false
|
||||
# Supported: sor, mg
|
||||
SOLVER ?= mg
|
||||
# Run in debug settings
|
||||
DEBUG ?= false
|
||||
|
||||
#Feature options
|
||||
OPTIONS += -DARRAY_ALIGNMENT=64
|
||||
#OPTIONS += -DVERBOSE
|
||||
#OPTIONS += -DDEBUG
|
||||
#OPTIONS += -DBOUNDCHECK
|
||||
#OPTIONS += -DVERBOSE_AFFINITY
|
||||
#OPTIONS += -DVERBOSE_DATASIZE
|
||||
#OPTIONS += -DVERBOSE_TIMER
|
||||
|
@ -26,8 +26,8 @@ p_init 0.0 # initial value for pressure
|
||||
|
||||
xlength 1.0 # domain size in x-direction
|
||||
ylength 1.0 # domain size in y-direction
|
||||
imax 100 # number of interior cells in x-direction
|
||||
jmax 100 # number of interior cells in y-direction
|
||||
imax 128 # number of interior cells in x-direction
|
||||
jmax 128 # number of interior cells in y-direction
|
||||
|
||||
# Time Data:
|
||||
# ---------
|
||||
@ -50,7 +50,7 @@ levels 5 # Multigrid levels
|
||||
# -----------------------
|
||||
|
||||
numberOfParticles 200
|
||||
startTime 100
|
||||
startTime 2.0
|
||||
injectTimePeriod 0.5
|
||||
writeTimePeriod 0.2
|
||||
|
||||
|
@ -2,16 +2,18 @@ CC = clang
|
||||
GCC = cc
|
||||
LINKER = $(CC)
|
||||
|
||||
ifeq ($(ENABLE_OPENMP),true)
|
||||
ifeq ($(strip $(ENABLE_OPENMP)),true)
|
||||
OPENMP = -fopenmp
|
||||
#OPENMP = -Xpreprocessor -fopenmp #required on Macos with homebrew libomp
|
||||
LIBS = # -lomp
|
||||
endif
|
||||
ifeq ($(strip $(DEBUG)),true)
|
||||
CFLAGS = -O0 -g -std=c17
|
||||
else
|
||||
CFLAGS = -O3 -std=c17 $(OPENMP)
|
||||
endif
|
||||
|
||||
VERSION = --version
|
||||
# CFLAGS = -O3 -std=c17 $(OPENMP)
|
||||
CFLAGS = -Ofast -std=c17
|
||||
#CFLAGS = -Ofast -fnt-store=aggressive -std=c99 $(OPENMP) #AMD CLANG
|
||||
LFLAGS = $(OPENMP) -lm
|
||||
DEFINES = -D_GNU_SOURCE# -DDEBUG
|
||||
DEFINES = -D_GNU_SOURCE
|
||||
INCLUDES =
|
||||
|
@ -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*/
|
@ -5,10 +5,13 @@
|
||||
* license that can be found in the LICENSE file.
|
||||
*/
|
||||
#include <errno.h>
|
||||
#include <stddef.h>
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
|
||||
void* allocate(int alignment, size_t bytesize)
|
||||
#include "allocate.h"
|
||||
|
||||
void* allocate(size_t alignment, size_t bytesize)
|
||||
{
|
||||
int errorCode;
|
||||
void* ptr;
|
||||
|
@ -8,6 +8,6 @@
|
||||
#define __ALLOCATE_H_
|
||||
#include <stdlib.h>
|
||||
|
||||
extern void* allocate(int alignment, size_t bytesize);
|
||||
extern void* allocate(size_t alignment, size_t bytesize);
|
||||
|
||||
#endif
|
||||
|
649
EnhancedSolver/2D-seq/src/discretization.c
Normal file
649
EnhancedSolver/2D-seq/src/discretization.c
Normal file
@ -0,0 +1,649 @@
|
||||
/*
|
||||
* Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
|
||||
* All rights reserved. This file is part of nusif-solver.
|
||||
* Use of this source code is governed by a MIT style
|
||||
* license that can be found in the LICENSE file.
|
||||
*/
|
||||
#include <float.h>
|
||||
#include <math.h>
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
|
||||
#include "allocate.h"
|
||||
#include "discretization.h"
|
||||
#include "grid.h"
|
||||
#include "parameter.h"
|
||||
#include "util.h"
|
||||
|
||||
#define S(i, j) s[(j) * (imax + 2) + (i)]
|
||||
|
||||
static double distance(double i, double j, double iCenter, double jCenter)
|
||||
{
|
||||
return sqrt(pow(iCenter - i, 2) + pow(jCenter - j, 2) * 1.0);
|
||||
}
|
||||
|
||||
void print(Discretization* d, double* grid)
|
||||
{
|
||||
int imax = d->grid.imax;
|
||||
int jmax = d->grid.jmax;
|
||||
|
||||
for (int j = 0; j < jmax + 2; j++) {
|
||||
printf("%02d: ", j);
|
||||
for (int i = 0; i < imax + 2; i++) {
|
||||
printf("%3.2f ", grid[j * (imax + 2) + i]);
|
||||
}
|
||||
printf("\n");
|
||||
}
|
||||
fflush(stdout);
|
||||
}
|
||||
|
||||
void printGrid(Discretization* d, int* grid)
|
||||
{
|
||||
int imax = d->grid.imax;
|
||||
int jmax = d->grid.jmax;
|
||||
|
||||
for (int j = 0; j < jmax + 2; j++) {
|
||||
printf("%02d: ", j);
|
||||
for (int i = 0; i < imax + 2; i++) {
|
||||
printf("%2d ", grid[j * (imax + 2) + i]);
|
||||
}
|
||||
printf("\n");
|
||||
}
|
||||
fflush(stdout);
|
||||
}
|
||||
|
||||
static void printConfig(Discretization* d)
|
||||
{
|
||||
printf("Parameters for #%s#\n", d->problem);
|
||||
printf("Boundary conditions Left:%d Right:%d Bottom:%d Top:%d\n",
|
||||
d->bcLeft,
|
||||
d->bcRight,
|
||||
d->bcBottom,
|
||||
d->bcTop);
|
||||
printf("\tReynolds number: %.2f\n", d->re);
|
||||
printf("\tGx Gy: %.2f %.2f\n", d->gx, d->gy);
|
||||
printf("Geometry data:\n");
|
||||
printf("\tDomain box size (x, y): %.2f, %.2f\n", d->grid.xlength, d->grid.ylength);
|
||||
printf("\tCells (x, y): %d, %d\n", d->grid.imax, d->grid.jmax);
|
||||
printf("Timestep parameters:\n");
|
||||
printf("\tDefault stepsize: %.2f, Final time %.2f\n", d->dt, d->te);
|
||||
printf("\tdt bound: %.6f\n", d->dtBound);
|
||||
printf("\tTau factor: %.2f\n", d->tau);
|
||||
printf("Iterative d parameters:\n");
|
||||
printf("\tgamma factor: %f\n", d->gamma);
|
||||
}
|
||||
|
||||
void initDiscretization(Discretization* d, Parameter* p)
|
||||
{
|
||||
d->problem = p->name;
|
||||
d->bcLeft = p->bcLeft;
|
||||
d->bcRight = p->bcRight;
|
||||
d->bcBottom = p->bcBottom;
|
||||
d->bcTop = p->bcTop;
|
||||
d->grid.imax = p->imax;
|
||||
d->grid.jmax = p->jmax;
|
||||
d->grid.xlength = p->xlength;
|
||||
d->grid.ylength = p->ylength;
|
||||
d->grid.dx = p->xlength / p->imax;
|
||||
d->grid.dy = p->ylength / p->jmax;
|
||||
d->re = p->re;
|
||||
d->gx = p->gx;
|
||||
d->gy = p->gy;
|
||||
d->dt = p->dt;
|
||||
d->te = p->te;
|
||||
d->tau = p->tau;
|
||||
d->gamma = p->gamma;
|
||||
|
||||
int imax = d->grid.imax;
|
||||
int jmax = d->grid.jmax;
|
||||
size_t size = (imax + 2) * (jmax + 2) * sizeof(double);
|
||||
d->u = allocate(64, size);
|
||||
d->v = allocate(64, size);
|
||||
d->grid.s = allocate(64, size);
|
||||
d->p = allocate(64, size);
|
||||
d->rhs = allocate(64, size);
|
||||
d->f = allocate(64, size);
|
||||
d->g = allocate(64, size);
|
||||
|
||||
for (int i = 0; i < (imax + 2) * (jmax + 2); i++) {
|
||||
|
||||
d->u[i] = p->u_init;
|
||||
d->v[i] = p->v_init;
|
||||
d->p[i] = p->p_init;
|
||||
d->rhs[i] = 0.0;
|
||||
d->f[i] = 0.0;
|
||||
d->g[i] = 0.0;
|
||||
d->grid.s[i] = FLUID;
|
||||
}
|
||||
|
||||
double dx = d->grid.dx;
|
||||
double dy = d->grid.dy;
|
||||
double invSqrSum = 1.0 / (dx * dx) + 1.0 / (dy * dy);
|
||||
d->dtBound = 0.5 * d->re * 1.0 / invSqrSum;
|
||||
|
||||
double xCenter = 0, yCenter = 0, radius = 0;
|
||||
double x1 = 0, x2 = 0, y1 = 0, y2 = 0;
|
||||
|
||||
int* s = d->grid.s;
|
||||
|
||||
switch (p->shape) {
|
||||
case NOSHAPE:
|
||||
break;
|
||||
case RECT:
|
||||
x1 = p->xCenter - p->xRectLength / 2;
|
||||
x2 = p->xCenter + p->xRectLength / 2;
|
||||
y1 = p->yCenter - p->yRectLength / 2;
|
||||
y2 = p->yCenter + p->yRectLength / 2;
|
||||
|
||||
for (int j = 1; j < jmax + 1; j++) {
|
||||
for (int i = 1; i < imax + 1; i++) {
|
||||
if ((x1 <= (i * dx)) && ((i * dx) <= x2) && (y1 <= (j * dy)) &&
|
||||
((j * dy) <= y2)) {
|
||||
S(i, j) = OBSTACLE;
|
||||
}
|
||||
}
|
||||
}
|
||||
break;
|
||||
case CIRCLE:
|
||||
xCenter = p->xCenter;
|
||||
yCenter = p->yCenter;
|
||||
radius = p->circleRadius;
|
||||
|
||||
for (int j = 1; j < jmax + 1; j++) {
|
||||
for (int i = 1; i < imax + 1; i++) {
|
||||
if (distance((i * dx), (j * dy), xCenter, yCenter) <= radius) {
|
||||
S(i, j) = OBSTACLE;
|
||||
}
|
||||
}
|
||||
}
|
||||
break;
|
||||
}
|
||||
|
||||
if (p->shape != NOSHAPE) {
|
||||
for (int j = 1; j < jmax + 1; j++) {
|
||||
for (int i = 1; i < imax + 1; i++) {
|
||||
|
||||
if (S(i, j - 1) == FLUID && S(i, j + 1) == OBSTACLE &&
|
||||
S(i, j) == OBSTACLE)
|
||||
S(i, j) = BOTTOM; // TOP
|
||||
if (S(i - 1, j) == FLUID && S(i + 1, j) == OBSTACLE &&
|
||||
S(i, j) == OBSTACLE)
|
||||
S(i, j) = LEFT; // LEFT
|
||||
if (S(i + 1, j) == FLUID && S(i - 1, j) == OBSTACLE &&
|
||||
S(i, j) == OBSTACLE)
|
||||
S(i, j) = RIGHT; // RIGHT
|
||||
if (S(i, j + 1) == FLUID && S(i, j - 1) == OBSTACLE &&
|
||||
S(i, j) == OBSTACLE)
|
||||
S(i, j) = TOP; // BOTTOM
|
||||
if (S(i - 1, j - 1) == FLUID && S(i, j - 1) == FLUID &&
|
||||
S(i - 1, j) == FLUID && S(i + 1, j + 1) == OBSTACLE &&
|
||||
(S(i, j) == OBSTACLE || S(i, j) == LEFT || S(i, j) == BOTTOM))
|
||||
S(i, j) = BOTTOMLEFT; // TOPLEFT
|
||||
if (S(i + 1, j - 1) == FLUID && S(i, j - 1) == FLUID &&
|
||||
S(i + 1, j) == FLUID && S(i - 1, j + 1) == OBSTACLE &&
|
||||
(S(i, j) == OBSTACLE || S(i, j) == RIGHT || S(i, j) == BOTTOM))
|
||||
S(i, j) = BOTTOMRIGHT; // TOPRIGHT
|
||||
if (S(i - 1, j + 1) == FLUID && S(i - 1, j) == FLUID &&
|
||||
S(i, j + 1) == FLUID && S(i + 1, j - 1) == OBSTACLE &&
|
||||
(S(i, j) == OBSTACLE || S(i, j) == LEFT || S(i, j) == TOP))
|
||||
S(i, j) = TOPLEFT; // BOTTOMLEFT
|
||||
if (S(i + 1, j + 1) == FLUID && S(i + 1, j) == FLUID &&
|
||||
S(i, j + 1) == FLUID && S(i - 1, j - 1) == OBSTACLE &&
|
||||
(S(i, j) == OBSTACLE || S(i, j) == RIGHT || S(i, j) == TOP))
|
||||
S(i, j) = TOPRIGHT; // BOTTOMRIGHT
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#ifdef VERBOSE
|
||||
printConfig(solver);
|
||||
#endif
|
||||
}
|
||||
|
||||
static double maxElement(Discretization* d, double* m)
|
||||
{
|
||||
int size = (d->grid.imax + 2) * (d->grid.jmax + 2);
|
||||
double maxval = DBL_MIN;
|
||||
|
||||
for (int i = 0; i < size; i++) {
|
||||
maxval = MAX(maxval, fabs(m[i]));
|
||||
}
|
||||
|
||||
return maxval;
|
||||
}
|
||||
|
||||
void computeRHS(Discretization* d)
|
||||
{
|
||||
int imax = d->grid.imax;
|
||||
int jmax = d->grid.jmax;
|
||||
double idx = 1.0 / d->grid.dx;
|
||||
double idy = 1.0 / d->grid.dy;
|
||||
double idt = 1.0 / d->dt;
|
||||
double* rhs = d->rhs;
|
||||
double* f = d->f;
|
||||
double* g = d->g;
|
||||
int* s = d->grid.s;
|
||||
|
||||
for (int j = 1; j < jmax + 1; j++) {
|
||||
for (int i = 1; i < imax + 1; i++) {
|
||||
RHS(i, j) = idt *
|
||||
((F(i, j) - F(i - 1, j)) * idx + (G(i, j) - G(i, j - 1)) * idy);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void normalizePressure(Discretization* d)
|
||||
{
|
||||
int size = (d->grid.imax + 2) * (d->grid.jmax + 2);
|
||||
double* p = d->p;
|
||||
double avgP = 0.0;
|
||||
|
||||
for (int i = 0; i < size; i++) {
|
||||
avgP += p[i];
|
||||
}
|
||||
avgP /= size;
|
||||
|
||||
for (int i = 0; i < size; i++) {
|
||||
p[i] = p[i] - avgP;
|
||||
}
|
||||
}
|
||||
|
||||
void computeTimestep(Discretization* d)
|
||||
{
|
||||
double dt = d->dtBound;
|
||||
double dx = d->grid.dx;
|
||||
double dy = d->grid.dy;
|
||||
double umax = maxElement(d, d->u);
|
||||
double vmax = maxElement(d, d->v);
|
||||
|
||||
if (umax > 0) {
|
||||
dt = (dt > dx / umax) ? dx / umax : dt;
|
||||
}
|
||||
if (vmax > 0) {
|
||||
dt = (dt > dy / vmax) ? dy / vmax : dt;
|
||||
}
|
||||
|
||||
d->dt = dt * d->tau;
|
||||
}
|
||||
|
||||
void setBoundaryConditions(Discretization* d)
|
||||
{
|
||||
int imax = d->grid.imax;
|
||||
int jmax = d->grid.jmax;
|
||||
double* u = d->u;
|
||||
double* v = d->v;
|
||||
|
||||
// Left boundary
|
||||
switch (d->bcLeft) {
|
||||
case NOSLIP:
|
||||
for (int j = 1; j < jmax + 1; j++) {
|
||||
U(0, j) = 0.0;
|
||||
V(0, j) = -V(1, j);
|
||||
}
|
||||
break;
|
||||
case SLIP:
|
||||
for (int j = 1; j < jmax + 1; j++) {
|
||||
U(0, j) = 0.0;
|
||||
V(0, j) = V(1, j);
|
||||
}
|
||||
break;
|
||||
case OUTFLOW:
|
||||
for (int j = 1; j < jmax + 1; j++) {
|
||||
U(0, j) = U(1, j);
|
||||
V(0, j) = V(1, j);
|
||||
}
|
||||
break;
|
||||
case PERIODIC:
|
||||
break;
|
||||
}
|
||||
|
||||
// Right boundary
|
||||
switch (d->bcRight) {
|
||||
case NOSLIP:
|
||||
for (int j = 1; j < jmax + 1; j++) {
|
||||
U(imax, j) = 0.0;
|
||||
V(imax + 1, j) = -V(imax, j);
|
||||
}
|
||||
break;
|
||||
case SLIP:
|
||||
for (int j = 1; j < jmax + 1; j++) {
|
||||
U(imax, j) = 0.0;
|
||||
V(imax + 1, j) = V(imax, j);
|
||||
}
|
||||
break;
|
||||
case OUTFLOW:
|
||||
for (int j = 1; j < jmax + 1; j++) {
|
||||
U(imax, j) = U(imax - 1, j);
|
||||
V(imax + 1, j) = V(imax, j);
|
||||
}
|
||||
break;
|
||||
case PERIODIC:
|
||||
break;
|
||||
}
|
||||
|
||||
// Bottom boundary
|
||||
switch (d->bcBottom) {
|
||||
case NOSLIP:
|
||||
for (int i = 1; i < imax + 1; i++) {
|
||||
V(i, 0) = 0.0;
|
||||
U(i, 0) = -U(i, 1);
|
||||
}
|
||||
break;
|
||||
case SLIP:
|
||||
for (int i = 1; i < imax + 1; i++) {
|
||||
V(i, 0) = 0.0;
|
||||
U(i, 0) = U(i, 1);
|
||||
}
|
||||
break;
|
||||
case OUTFLOW:
|
||||
for (int i = 1; i < imax + 1; i++) {
|
||||
U(i, 0) = U(i, 1);
|
||||
V(i, 0) = V(i, 1);
|
||||
}
|
||||
break;
|
||||
case PERIODIC:
|
||||
break;
|
||||
}
|
||||
|
||||
// Top boundary
|
||||
switch (d->bcTop) {
|
||||
case NOSLIP:
|
||||
for (int i = 1; i < imax + 1; i++) {
|
||||
V(i, jmax) = 0.0;
|
||||
U(i, jmax + 1) = -U(i, jmax);
|
||||
}
|
||||
break;
|
||||
case SLIP:
|
||||
for (int i = 1; i < imax + 1; i++) {
|
||||
V(i, jmax) = 0.0;
|
||||
U(i, jmax + 1) = U(i, jmax);
|
||||
}
|
||||
break;
|
||||
case OUTFLOW:
|
||||
for (int i = 1; i < imax + 1; i++) {
|
||||
U(i, jmax + 1) = U(i, jmax);
|
||||
V(i, jmax) = V(i, jmax - 1);
|
||||
}
|
||||
break;
|
||||
case PERIODIC:
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
void setSpecialBoundaryCondition(Discretization* d)
|
||||
{
|
||||
int imax = d->grid.imax;
|
||||
int jmax = d->grid.jmax;
|
||||
double mDy = d->grid.dy;
|
||||
double* u = d->u;
|
||||
int* s = d->grid.s;
|
||||
|
||||
if (strcmp(d->problem, "dcavity") == 0) {
|
||||
for (int i = 1; i < imax; i++) {
|
||||
U(i, jmax + 1) = 2.0 - U(i, jmax);
|
||||
}
|
||||
} else if (strcmp(d->problem, "canal") == 0) {
|
||||
double ylength = d->grid.ylength;
|
||||
double y;
|
||||
|
||||
for (int j = 1; j < jmax + 1; j++) {
|
||||
y = mDy * (j - 0.5);
|
||||
U(0, j) = y * (ylength - y) * 4.0 / (ylength * ylength);
|
||||
}
|
||||
} else if (strcmp(d->problem, "backstep") == 0) {
|
||||
for (int j = 1; j < jmax + 1; j++) {
|
||||
if (S(0, j) == FLUID) U(0, j) = 1.0;
|
||||
}
|
||||
} else if (strcmp(d->problem, "karman") == 0) {
|
||||
for (int j = 1; j < jmax + 1; j++) {
|
||||
U(0, j) = 1.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void setObjectBoundaryCondition(Discretization* d)
|
||||
{
|
||||
int imax = d->grid.imax;
|
||||
int jmax = d->grid.jmax;
|
||||
double* u = d->u;
|
||||
double* v = d->v;
|
||||
int* s = d->grid.s;
|
||||
|
||||
for (int j = 1; j < jmax + 1; j++) {
|
||||
for (int i = 1; i < imax + 1; i++) {
|
||||
switch (S(i, j)) {
|
||||
case TOP:
|
||||
U(i, j) = -U(i, j + 1);
|
||||
U(i - 1, j) = -U(i - 1, j + 1);
|
||||
V(i, j) = 0.0;
|
||||
break;
|
||||
case BOTTOM:
|
||||
U(i, j) = -U(i, j - 1);
|
||||
U(i - 1, j) = -U(i - 1, j - 1);
|
||||
V(i, j) = 0.0;
|
||||
break;
|
||||
case LEFT:
|
||||
U(i - 1, j) = 0.0;
|
||||
V(i, j) = -V(i - 1, j);
|
||||
V(i, j - 1) = -V(i - 1, j - 1);
|
||||
break;
|
||||
case RIGHT:
|
||||
U(i, j) = 0.0;
|
||||
V(i, j) = -V(i + 1, j);
|
||||
V(i, j - 1) = -V(i + 1, j - 1);
|
||||
break;
|
||||
case TOPLEFT:
|
||||
U(i, j) = -U(i, j + 1);
|
||||
U(i - 1, j) = 0.0;
|
||||
V(i, j) = 0.0;
|
||||
V(i, j - 1) = -V(i - 1, j - 1);
|
||||
break;
|
||||
case TOPRIGHT:
|
||||
U(i, j) = 0.0;
|
||||
U(i - 1, j) = -U(i - 1, j + 1);
|
||||
V(i, j) = 0.0;
|
||||
V(i, j - 1) = -V(i + 1, j - 1);
|
||||
break;
|
||||
case BOTTOMLEFT:
|
||||
U(i, j) = -U(i, j - 1);
|
||||
U(i - 1, j) = 0.0;
|
||||
V(i, j) = -V(i - 1, j);
|
||||
V(i, j - 1) = 0.0;
|
||||
break;
|
||||
case BOTTOMRIGHT:
|
||||
U(i, j) = 0.0;
|
||||
U(i - 1, j) = -U(i - 1, j - 1);
|
||||
V(i, j) = -V(i, j + 1);
|
||||
V(i, j - 1) = 0.0;
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void computeFG(Discretization* d)
|
||||
{
|
||||
double* u = d->u;
|
||||
double* v = d->v;
|
||||
double* f = d->f;
|
||||
double* g = d->g;
|
||||
int* s = d->grid.s;
|
||||
int imax = d->grid.imax;
|
||||
int jmax = d->grid.jmax;
|
||||
double gx = d->gx;
|
||||
double gy = d->gy;
|
||||
double gamma = d->gamma;
|
||||
double dt = d->dt;
|
||||
double inverseRe = 1.0 / d->re;
|
||||
double inverseDx = 1.0 / d->grid.dx;
|
||||
double inverseDy = 1.0 / d->grid.dy;
|
||||
double du2dx, dv2dy, duvdx, duvdy;
|
||||
double du2dx2, du2dy2, dv2dx2, dv2dy2;
|
||||
|
||||
for (int j = 1; j < jmax + 1; j++) {
|
||||
for (int i = 1; i < imax + 1; i++) {
|
||||
if (S(i, j) == FLUID) {
|
||||
du2dx = inverseDx * 0.25 *
|
||||
((U(i, j) + U(i + 1, j)) * (U(i, j) + U(i + 1, j)) -
|
||||
(U(i, j) + U(i - 1, j)) * (U(i, j) + U(i - 1, j))) +
|
||||
gamma * inverseDx * 0.25 *
|
||||
(fabs(U(i, j) + U(i + 1, j)) * (U(i, j) - U(i + 1, j)) +
|
||||
fabs(U(i, j) + U(i - 1, j)) * (U(i, j) - U(i - 1, j)));
|
||||
|
||||
duvdy = inverseDy * 0.25 *
|
||||
((V(i, j) + V(i + 1, j)) * (U(i, j) + U(i, j + 1)) -
|
||||
(V(i, j - 1) + V(i + 1, j - 1)) *
|
||||
(U(i, j) + U(i, j - 1))) +
|
||||
gamma * inverseDy * 0.25 *
|
||||
(fabs(V(i, j) + V(i + 1, j)) * (U(i, j) - U(i, j + 1)) +
|
||||
fabs(V(i, j - 1) + V(i + 1, j - 1)) *
|
||||
(U(i, j) - U(i, j - 1)));
|
||||
|
||||
du2dx2 = inverseDx * inverseDx *
|
||||
(U(i + 1, j) - 2.0 * U(i, j) + U(i - 1, j));
|
||||
du2dy2 = inverseDy * inverseDy *
|
||||
(U(i, j + 1) - 2.0 * U(i, j) + U(i, j - 1));
|
||||
F(i, j) = U(i, j) +
|
||||
dt * (inverseRe * (du2dx2 + du2dy2) - du2dx - duvdy + gx);
|
||||
|
||||
duvdx = inverseDx * 0.25 *
|
||||
((U(i, j) + U(i, j + 1)) * (V(i, j) + V(i + 1, j)) -
|
||||
(U(i - 1, j) + U(i - 1, j + 1)) *
|
||||
(V(i, j) + V(i - 1, j))) +
|
||||
gamma * inverseDx * 0.25 *
|
||||
(fabs(U(i, j) + U(i, j + 1)) * (V(i, j) - V(i + 1, j)) +
|
||||
fabs(U(i - 1, j) + U(i - 1, j + 1)) *
|
||||
(V(i, j) - V(i - 1, j)));
|
||||
|
||||
dv2dy = inverseDy * 0.25 *
|
||||
((V(i, j) + V(i, j + 1)) * (V(i, j) + V(i, j + 1)) -
|
||||
(V(i, j) + V(i, j - 1)) * (V(i, j) + V(i, j - 1))) +
|
||||
gamma * inverseDy * 0.25 *
|
||||
(fabs(V(i, j) + V(i, j + 1)) * (V(i, j) - V(i, j + 1)) +
|
||||
fabs(V(i, j) + V(i, j - 1)) * (V(i, j) - V(i, j - 1)));
|
||||
|
||||
dv2dx2 = inverseDx * inverseDx *
|
||||
(V(i + 1, j) - 2.0 * V(i, j) + V(i - 1, j));
|
||||
dv2dy2 = inverseDy * inverseDy *
|
||||
(V(i, j + 1) - 2.0 * V(i, j) + V(i, j - 1));
|
||||
G(i, j) = V(i, j) +
|
||||
dt * (inverseRe * (dv2dx2 + dv2dy2) - duvdx - dv2dy + gy);
|
||||
} else {
|
||||
switch (S(i, j)) {
|
||||
case TOP:
|
||||
G(i, j) = V(i, j);
|
||||
break;
|
||||
case BOTTOM:
|
||||
G(i, j - 1) = V(i, j - 1);
|
||||
break;
|
||||
case LEFT:
|
||||
F(i - 1, j) = U(i - 1, j);
|
||||
break;
|
||||
case RIGHT:
|
||||
F(i, j) = U(i, j);
|
||||
break;
|
||||
case TOPLEFT:
|
||||
F(i - 1, j) = U(i - 1, j);
|
||||
G(i, j) = V(i, j);
|
||||
break;
|
||||
case TOPRIGHT:
|
||||
F(i, j) = U(i, j);
|
||||
G(i, j) = V(i, j);
|
||||
break;
|
||||
case BOTTOMLEFT:
|
||||
F(i - 1, j) = U(i - 1, j);
|
||||
G(i, j - 1) = V(i, j - 1);
|
||||
break;
|
||||
case BOTTOMRIGHT:
|
||||
F(i, j) = U(i, j);
|
||||
G(i, j - 1) = V(i, j - 1);
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------- boundary of F --------------------------- */
|
||||
for (int j = 1; j < jmax + 1; j++) {
|
||||
F(0, j) = U(0, j);
|
||||
F(imax, j) = U(imax, j);
|
||||
}
|
||||
|
||||
/* ---------------------- boundary of G --------------------------- */
|
||||
for (int i = 1; i < imax + 1; i++) {
|
||||
G(i, 0) = V(i, 0);
|
||||
G(i, jmax) = V(i, jmax);
|
||||
}
|
||||
}
|
||||
|
||||
void adaptUV(Discretization* d)
|
||||
{
|
||||
int imax = d->grid.imax;
|
||||
int jmax = d->grid.jmax;
|
||||
double* p = d->p;
|
||||
double* u = d->u;
|
||||
double* v = d->v;
|
||||
double* f = d->f;
|
||||
double* g = d->g;
|
||||
double factorX = d->dt / d->grid.dx;
|
||||
double factorY = d->dt / d->grid.dy;
|
||||
|
||||
for (int j = 1; j < jmax + 1; j++) {
|
||||
for (int i = 1; i < imax + 1; i++) {
|
||||
U(i, j) = F(i, j) - (P(i + 1, j) - P(i, j)) * factorX;
|
||||
V(i, j) = G(i, j) - (P(i, j + 1) - P(i, j)) * factorY;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void writeResult(Discretization* d)
|
||||
{
|
||||
int imax = d->grid.imax;
|
||||
int jmax = d->grid.jmax;
|
||||
double dx = d->grid.dx;
|
||||
double dy = d->grid.dy;
|
||||
double* p = d->p;
|
||||
double* u = d->u;
|
||||
double* v = d->v;
|
||||
double x = 0.0, y = 0.0;
|
||||
|
||||
FILE* fp;
|
||||
fp = fopen("pressure.dat", "w");
|
||||
|
||||
if (fp == NULL) {
|
||||
printf("Error!\n");
|
||||
exit(EXIT_FAILURE);
|
||||
}
|
||||
|
||||
for (int j = 1; j < jmax + 1; j++) {
|
||||
y = (double)(j - 0.5) * dy;
|
||||
for (int i = 1; i < imax + 1; i++) {
|
||||
x = (double)(i - 0.5) * dx;
|
||||
fprintf(fp, "%.2f %.2f %f\n", x, y, P(i, j));
|
||||
}
|
||||
fprintf(fp, "\n");
|
||||
}
|
||||
|
||||
fclose(fp);
|
||||
|
||||
fp = fopen("velocity.dat", "w");
|
||||
|
||||
if (fp == NULL) {
|
||||
printf("Error!\n");
|
||||
exit(EXIT_FAILURE);
|
||||
}
|
||||
|
||||
for (int j = 1; j < jmax + 1; j++) {
|
||||
y = dy * (j - 0.5);
|
||||
for (int i = 1; i < imax + 1; i++) {
|
||||
x = dx * (i - 0.5);
|
||||
double velU = (U(i, j) + U(i - 1, j)) / 2.0;
|
||||
double velV = (V(i, j) + V(i, j - 1)) / 2.0;
|
||||
double len = sqrt((velU * velU) + (velV * velV));
|
||||
fprintf(fp, "%.2f %.2f %f %f %f\n", x, y, velU, velV, len);
|
||||
}
|
||||
}
|
||||
|
||||
fclose(fp);
|
||||
}
|
44
EnhancedSolver/2D-seq/src/discretization.h
Normal file
44
EnhancedSolver/2D-seq/src/discretization.h
Normal file
@ -0,0 +1,44 @@
|
||||
/*
|
||||
* Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
|
||||
* All rights reserved. This file is part of nusif-solver.
|
||||
* Use of this source code is governed by a MIT style
|
||||
* license that can be found in the LICENSE file.
|
||||
*/
|
||||
#ifndef __DISCRETIZATION_H_
|
||||
#define __DISCRETIZATION_H_
|
||||
#include "grid.h"
|
||||
#include "parameter.h"
|
||||
|
||||
enum BC { NOSLIP = 1, SLIP, OUTFLOW, PERIODIC };
|
||||
|
||||
typedef struct {
|
||||
/* geometry and grid information */
|
||||
Grid grid;
|
||||
/* arrays */
|
||||
double *p, *rhs;
|
||||
double *f, *g;
|
||||
double *u, *v;
|
||||
/* parameters */
|
||||
double rho;
|
||||
double re, tau, gamma;
|
||||
double gx, gy;
|
||||
/* time stepping */
|
||||
double dt, te;
|
||||
double dtBound;
|
||||
char* problem;
|
||||
int bcLeft, bcRight, bcBottom, bcTop;
|
||||
} Discretization;
|
||||
|
||||
extern void initDiscretization(Discretization*, Parameter*);
|
||||
extern void computeRHS(Discretization*);
|
||||
extern void normalizePressure(Discretization*);
|
||||
extern void computeTimestep(Discretization*);
|
||||
extern void setBoundaryConditions(Discretization*);
|
||||
extern void setSpecialBoundaryCondition(Discretization*);
|
||||
extern void setObjectBoundaryCondition(Discretization*);
|
||||
extern void computeFG(Discretization*);
|
||||
extern void adaptUV(Discretization*);
|
||||
extern void writeResult(Discretization*);
|
||||
extern void print(Discretization*, double*);
|
||||
extern void printGrid(Discretization*, int*);
|
||||
#endif
|
@ -7,10 +7,32 @@
|
||||
#ifndef __GRID_H_
|
||||
#define __GRID_H_
|
||||
|
||||
#define S(i, j) s[(j) * (imax + 2) + (i)]
|
||||
|
||||
enum OBJECTBOUNDARY {
|
||||
FLUID = 0,
|
||||
TOP,
|
||||
BOTTOM,
|
||||
LEFT,
|
||||
RIGHT,
|
||||
TOPLEFT,
|
||||
BOTTOMLEFT,
|
||||
TOPRIGHT,
|
||||
BOTTOMRIGHT,
|
||||
OBSTACLE
|
||||
};
|
||||
|
||||
enum SHAPE { NOSHAPE = 0, RECT, CIRCLE };
|
||||
|
||||
typedef struct {
|
||||
double dx, dy;
|
||||
int imax, jmax;
|
||||
double xlength, ylength;
|
||||
int* s;
|
||||
} Grid;
|
||||
|
||||
static inline int gridIsFluid(Grid* g, int i, int j)
|
||||
{
|
||||
return g->s[j * (g->imax + 2) + i];
|
||||
}
|
||||
#endif // __GRID_H_
|
||||
|
@ -8,6 +8,7 @@
|
||||
#include <stdlib.h>
|
||||
#include <unistd.h>
|
||||
|
||||
#include "discretization.h"
|
||||
#include "parameter.h"
|
||||
#include "particletracing.h"
|
||||
#include "progress.h"
|
||||
@ -16,52 +17,51 @@
|
||||
|
||||
int main(int argc, char** argv)
|
||||
{
|
||||
double S, E;
|
||||
Parameter params;
|
||||
Solver solver;
|
||||
double timeStart, timeStop;
|
||||
Parameter p;
|
||||
Discretization d;
|
||||
Solver s;
|
||||
ParticleTracer particletracer;
|
||||
initParameter(¶ms);
|
||||
|
||||
initParameter(&p);
|
||||
|
||||
if (argc != 2) {
|
||||
printf("Usage: %s <configFile>\n", argv[0]);
|
||||
exit(EXIT_SUCCESS);
|
||||
}
|
||||
|
||||
readParameter(¶ms, argv[1]);
|
||||
printParameter(¶ms);
|
||||
initSolver(&solver, ¶ms);
|
||||
printf("initsolver done\n");
|
||||
|
||||
initParticleTracer(&particletracer, &solver.grid, ¶ms);
|
||||
readParameter(&p, argv[1]);
|
||||
printParameter(&p);
|
||||
initDiscretization(&d, &p);
|
||||
initSolver(&s, &d, &p);
|
||||
initParticleTracer(&particletracer, &d.grid, &p);
|
||||
printParticleTracerParameters(&particletracer);
|
||||
|
||||
#ifndef VERBOSE
|
||||
initProgress(solver.te);
|
||||
initProgress(d.te);
|
||||
#endif
|
||||
|
||||
double tau = solver.tau;
|
||||
double te = solver.te;
|
||||
double tau = d.tau;
|
||||
double te = d.te;
|
||||
double t = 0.0;
|
||||
int nt = 0;
|
||||
|
||||
S = getTimeStamp();
|
||||
timeStart = getTimeStamp();
|
||||
|
||||
while (t <= te) {
|
||||
if (tau > 0.0) computeTimestep(&solver);
|
||||
setBoundaryConditions(&solver);
|
||||
setSpecialBoundaryCondition(&solver);
|
||||
setObjectBoundaryCondition(&solver);
|
||||
computeFG(&solver);
|
||||
computeRHS(&solver);
|
||||
if (nt % 100 == 0) normalizePressure(&solver);
|
||||
multiGrid(&solver);
|
||||
if (tau > 0.0) computeTimestep(&d);
|
||||
setBoundaryConditions(&d);
|
||||
setSpecialBoundaryCondition(&d);
|
||||
setObjectBoundaryCondition(&d);
|
||||
|
||||
adaptUV(&solver);
|
||||
computeFG(&d);
|
||||
computeRHS(&d);
|
||||
if (nt % 100 == 0) normalizePressure(&d);
|
||||
solve(&s, d.p, d.rhs);
|
||||
adaptUV(&d);
|
||||
trace(&particletracer, d.u, d.v, d.dt, t);
|
||||
|
||||
/* Added function for particle tracing. Will inject and advance particles as per
|
||||
* timePeriod */
|
||||
trace(&particletracer, solver.u, solver.v, solver.s, t);
|
||||
|
||||
t += solver.dt;
|
||||
t += d.dt;
|
||||
nt++;
|
||||
|
||||
#ifdef VERBOSE
|
||||
@ -70,14 +70,12 @@ int main(int argc, char** argv)
|
||||
printProgress(t);
|
||||
#endif
|
||||
}
|
||||
printf("Total particles : %d\n", particletracer.totalParticles);
|
||||
|
||||
E = getTimeStamp();
|
||||
timeStop = getTimeStamp();
|
||||
|
||||
stopProgress();
|
||||
|
||||
freeParticles(&particletracer);
|
||||
|
||||
printf("Solution took %.2fs\n", E - S);
|
||||
writeResult(&solver);
|
||||
printf("Solution took %.2fs\n", timeStop - timeStart);
|
||||
writeResult(&d);
|
||||
return EXIT_SUCCESS;
|
||||
}
|
||||
|
@ -8,14 +8,13 @@
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
|
||||
#include "grid.h"
|
||||
#include "vtkWriter.h"
|
||||
|
||||
#define U(i, j) u[(j) * (imax + 2) + (i)]
|
||||
#define V(i, j) v[(j) * (imax + 2) + (i)]
|
||||
#define S(i, j) s[(j) * (imax + 2) + (i)]
|
||||
|
||||
static int ts = 0;
|
||||
|
||||
void printParticles(ParticleTracer* p)
|
||||
{
|
||||
for (int i = 0; i < p->totalParticles; ++i) {
|
||||
@ -25,29 +24,28 @@ void printParticles(ParticleTracer* p)
|
||||
p->particlePool[i].flag);
|
||||
}
|
||||
}
|
||||
void injectParticles(ParticleTracer* p)
|
||||
|
||||
static void injectParticles(ParticleTracer* p)
|
||||
{
|
||||
for (int i = 0; i < p->numberOfParticles; ++i) {
|
||||
if (p->totalParticles + p->numParticlesInLine > p->numAllocatedParticles) {
|
||||
return;
|
||||
}
|
||||
for (int i = 0; i < p->numParticlesInLine; ++i) {
|
||||
p->particlePool[p->pointer].x = p->linSpaceLine[i].x;
|
||||
p->particlePool[p->pointer].y = p->linSpaceLine[i].y;
|
||||
p->particlePool[p->pointer].flag = true;
|
||||
++(p->pointer);
|
||||
++(p->totalParticles);
|
||||
p->pointer++;
|
||||
p->totalParticles++;
|
||||
}
|
||||
}
|
||||
|
||||
void advanceParticles(ParticleTracer* p,
|
||||
double* restrict u,
|
||||
double* restrict v,
|
||||
int* restrict s,
|
||||
double time)
|
||||
static void advanceParticles(
|
||||
ParticleTracer* p, double* restrict u, double* restrict v, double dt)
|
||||
{
|
||||
int imax = p->grid->imax;
|
||||
int jmax = p->grid->jmax;
|
||||
|
||||
double dx = p->grid->dx;
|
||||
double dy = p->grid->dy;
|
||||
|
||||
double xlength = p->grid->xlength;
|
||||
double ylength = p->grid->ylength;
|
||||
|
||||
@ -64,14 +62,14 @@ void advanceParticles(ParticleTracer* p,
|
||||
double x2 = (double)iCoord * dx;
|
||||
double y2 = ((double)jCoord - 0.5) * dy;
|
||||
|
||||
double u_n = (1.0 / (dx * dy)) *
|
||||
double intU = (1.0 / (dx * dy)) *
|
||||
((x2 - x) * (y2 - y) * U(iCoord - 1, jCoord - 1) +
|
||||
(x - x1) * (y2 - y) * U(iCoord, jCoord - 1) +
|
||||
(x2 - x) * (y - y1) * U(iCoord - 1, jCoord) +
|
||||
(x - x1) * (y - y1) * U(iCoord, jCoord));
|
||||
|
||||
double new_x = x + p->dt * u_n;
|
||||
p->particlePool[i].x = new_x;
|
||||
double newX = x + dt * intU;
|
||||
p->particlePool[i].x = newX;
|
||||
|
||||
iCoord = (int)((x + 0.5 * dx) / dx) + 1;
|
||||
jCoord = (int)(y / dy) + 1;
|
||||
@ -81,50 +79,55 @@ void advanceParticles(ParticleTracer* p,
|
||||
x2 = ((double)iCoord - 0.5) * dx;
|
||||
y2 = (double)jCoord * dy;
|
||||
|
||||
double v_n = (1.0 / (dx * dy)) *
|
||||
double intV = (1.0 / (dx * dy)) *
|
||||
((x2 - x) * (y2 - y) * V(iCoord - 1, jCoord - 1) +
|
||||
(x - x1) * (y2 - y) * V(iCoord, jCoord - 1) +
|
||||
(x2 - x) * (y - y1) * V(iCoord - 1, jCoord) +
|
||||
(x - x1) * (y - y1) * V(iCoord, jCoord));
|
||||
|
||||
double new_y = y + p->dt * v_n;
|
||||
p->particlePool[i].y = new_y;
|
||||
double newY = y + dt * intV;
|
||||
p->particlePool[i].y = newY;
|
||||
|
||||
// printf("\tOld X : %.2f, New X : %.2f, iCoord : %d\n\tOld Y : %.2f, New Y :
|
||||
// %.2f, jCoord : %d\n\n", x, new_x, iCoord, y, new_y, jCoord);
|
||||
// printf("\tU(iCoord - 1, jCoord - 1) : %.2f, U(iCoord, jCoord - 1) : %.2f,
|
||||
// U(iCoord - 1, jCoord) : %.2f, U(iCoord, jCoord) : %.2f\n", U(iCoord - 1,
|
||||
// jCoord - 1), U(iCoord, jCoord - 1), U(iCoord - 1, jCoord), U(iCoord,
|
||||
// jCoord)); printf("\tV(iCoord - 1, jCoord - 1) : %.2f, V(iCoord, jCoord - 1)
|
||||
// : %.2f, V(iCoord - 1, jCoord) : %.2f, V(iCoord, jCoord) : %.2f\n\n",
|
||||
// V(iCoord - 1, jCoord - 1), V(iCoord, jCoord - 1), V(iCoord - 1, jCoord),
|
||||
// V(iCoord, jCoord)); printf("\t U N : %.2f, V N : %.2f\n\n", u_n, v_n);
|
||||
// printf("\t j-1 * (imax + 2) + i-1 = %d with element from U : %.2f", (jCoord
|
||||
// - 1) * (200 + 2) + (iCoord - 1), u[(jCoord - 1) * (imax + 2) + (iCoord -
|
||||
// 1)]); printf("\nimax : %d, jmax : %d\n", imax, jmax);
|
||||
|
||||
if (((new_x < 0.0) || (new_x > xlength) || (new_y < 0.0) ||
|
||||
(new_y > ylength))) {
|
||||
if (((newX < 0.0) || (newX > xlength) || (newY < 0.0) || (newY > ylength))) {
|
||||
p->particlePool[i].flag = false;
|
||||
p->removedParticles++;
|
||||
}
|
||||
int i_new = new_x / dx, j_new = new_y / dy;
|
||||
if (S(i_new, j_new) != NONE) {
|
||||
|
||||
int newI = newX / dx, newJ = newY / dy;
|
||||
|
||||
if (!gridIsFluid(p->grid, newI, newJ)) {
|
||||
p->particlePool[i].flag = false;
|
||||
p->removedParticles++;
|
||||
printf("Forbidden movement of particle into obstacle!\n");
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void freeParticles(ParticleTracer* p)
|
||||
static void compress(ParticleTracer* p)
|
||||
{
|
||||
if (p->particlePool != NULL) {
|
||||
free(p->particlePool);
|
||||
free(p->linSpaceLine);
|
||||
Particle* memPool = p->particlePool;
|
||||
Particle tempPool[p->totalParticles];
|
||||
int totalParticles = 0;
|
||||
|
||||
for (int i = 0; i < p->totalParticles; i++) {
|
||||
if (memPool[i].flag == 1) {
|
||||
tempPool[totalParticles].x = memPool[i].x;
|
||||
tempPool[totalParticles].y = memPool[i].y;
|
||||
tempPool[totalParticles].flag = memPool[i].flag;
|
||||
totalParticles++;
|
||||
}
|
||||
}
|
||||
|
||||
p->totalParticles = totalParticles;
|
||||
p->removedParticles = 0;
|
||||
p->pointer = totalParticles + 1;
|
||||
memcpy(p->particlePool, tempPool, totalParticles * sizeof(Particle));
|
||||
}
|
||||
|
||||
void writeParticles(ParticleTracer* p)
|
||||
{
|
||||
static int ts = 0;
|
||||
VtkOptions opts = { .particletracer = p };
|
||||
|
||||
char filename[50];
|
||||
@ -154,44 +157,43 @@ void writeParticles(ParticleTracer* p)
|
||||
++ts;
|
||||
}
|
||||
|
||||
void initParticleTracer(ParticleTracer* p, Grid* grid, Parameter* params)
|
||||
void initParticleTracer(ParticleTracer* pt, Grid* g, Parameter* p)
|
||||
{
|
||||
p->numberOfParticles = params->numberOfParticles;
|
||||
p->startTime = params->startTime;
|
||||
p->injectTimePeriod = params->injectTimePeriod;
|
||||
p->writeTimePeriod = params->writeTimePeriod;
|
||||
pt->numParticlesInLine = p->numberOfParticles;
|
||||
pt->startTime = p->startTime;
|
||||
pt->injectTimePeriod = p->injectTimePeriod;
|
||||
pt->writeTimePeriod = p->writeTimePeriod;
|
||||
pt->grid = g;
|
||||
|
||||
p->dt = params->dt;
|
||||
p->grid = grid;
|
||||
pt->x1 = p->x1;
|
||||
pt->y1 = p->y1;
|
||||
pt->x2 = p->x2;
|
||||
pt->y2 = p->y2;
|
||||
|
||||
p->x1 = params->x1;
|
||||
p->y1 = params->y1;
|
||||
p->x2 = params->x2;
|
||||
p->y2 = params->y2;
|
||||
pt->lastInjectTime = p->startTime;
|
||||
pt->lastWriteTime = p->startTime;
|
||||
|
||||
p->lastInjectTime = params->startTime;
|
||||
p->lastUpdateTime = params->startTime;
|
||||
p->lastWriteTime = params->startTime;
|
||||
pt->pointer = 0;
|
||||
pt->removedParticles = 0;
|
||||
pt->totalParticles = 0;
|
||||
|
||||
p->pointer = 0;
|
||||
p->totalParticles = 0;
|
||||
if (p->te > p->startTime) {
|
||||
pt->numAllocatedParticles = ((p->te - p->startTime) / p->injectTimePeriod) *
|
||||
p->numberOfParticles;
|
||||
pt->numAllocatedParticles += (2 * p->numberOfParticles);
|
||||
|
||||
if (params->te > params->startTime) {
|
||||
p->estimatedNumParticles = ((params->te - params->startTime) + 2) *
|
||||
params->numberOfParticles;
|
||||
pt->particlePool = malloc(sizeof(Particle) * pt->numAllocatedParticles);
|
||||
pt->linSpaceLine = malloc(sizeof(Particle) * pt->numParticlesInLine);
|
||||
|
||||
p->particlePool = malloc(sizeof(Particle) * p->estimatedNumParticles);
|
||||
p->linSpaceLine = malloc(sizeof(Particle) * p->numberOfParticles);
|
||||
|
||||
for (int i = 0; i < p->numberOfParticles; ++i) {
|
||||
double spacing = (double)i / (double)(p->numberOfParticles - 1);
|
||||
p->linSpaceLine[i].x = spacing * p->x1 + (1.0 - spacing) * p->x2;
|
||||
p->linSpaceLine[i].y = spacing * p->y1 + (1.0 - spacing) * p->y2;
|
||||
p->linSpaceLine[i].flag = true;
|
||||
for (int i = 0; i < pt->numParticlesInLine; ++i) {
|
||||
double spacing = (double)i / (double)(pt->numParticlesInLine - 1);
|
||||
pt->linSpaceLine[i].x = spacing * pt->x1 + (1.0 - spacing) * pt->x2;
|
||||
pt->linSpaceLine[i].y = spacing * pt->y1 + (1.0 - spacing) * pt->y2;
|
||||
pt->linSpaceLine[i].flag = true;
|
||||
}
|
||||
} else {
|
||||
p->particlePool = NULL;
|
||||
p->linSpaceLine = NULL;
|
||||
pt->particlePool = NULL;
|
||||
pt->linSpaceLine = NULL;
|
||||
}
|
||||
}
|
||||
|
||||
@ -199,7 +201,7 @@ void printParticleTracerParameters(ParticleTracer* p)
|
||||
{
|
||||
printf("Particle Tracing data:\n");
|
||||
printf("\tNumber of particles : %d being injected for every period of %.2f\n",
|
||||
p->numberOfParticles,
|
||||
p->numParticlesInLine,
|
||||
p->injectTimePeriod);
|
||||
printf("\tstartTime : %.2f\n", p->startTime);
|
||||
printf("\t(Line along which the particles are to be injected) \n\tx1 : %.2f, y1 : "
|
||||
@ -209,43 +211,33 @@ void printParticleTracerParameters(ParticleTracer* p)
|
||||
p->x2,
|
||||
p->y2);
|
||||
printf("\tPointer : %d, TotalParticles : %d\n", p->pointer, p->totalParticles);
|
||||
printf("\tdt : %.2f, dx : %.2f, dy : %.2f\n", p->dt, p->grid->dx, p->grid->dy);
|
||||
}
|
||||
|
||||
void trace(ParticleTracer* p, double* u, double* v, int* s, double time)
|
||||
void trace(ParticleTracer* p, double* u, double* v, double dt, double time)
|
||||
{
|
||||
if (time >= p->startTime) {
|
||||
// printParticles(particletracer);
|
||||
if ((time - p->lastInjectTime) >= p->injectTimePeriod) {
|
||||
injectParticles(p);
|
||||
p->lastInjectTime = time;
|
||||
}
|
||||
|
||||
if ((time - p->lastWriteTime) >= p->writeTimePeriod) {
|
||||
writeParticles(p);
|
||||
p->lastWriteTime = time;
|
||||
}
|
||||
advanceParticles(p, u, v, s, time);
|
||||
|
||||
advanceParticles(p, u, v, dt);
|
||||
|
||||
if (p->removedParticles > (p->totalParticles * 0.2)) {
|
||||
compress(p);
|
||||
p->lastUpdateTime = time;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void compress(ParticleTracer* p)
|
||||
void freeParticles(ParticleTracer* p)
|
||||
{
|
||||
Particle* memPool = p->particlePool;
|
||||
Particle tempPool[p->totalParticles];
|
||||
int totalParticles = 0;
|
||||
|
||||
for (int i = 0; i < p->totalParticles; ++i) {
|
||||
if (memPool[i].flag == 1) {
|
||||
tempPool[totalParticles].x = memPool[i].x;
|
||||
tempPool[totalParticles].y = memPool[i].y;
|
||||
tempPool[totalParticles].flag = memPool[i].flag;
|
||||
++totalParticles;
|
||||
if (p->particlePool != NULL) {
|
||||
free(p->particlePool);
|
||||
free(p->linSpaceLine);
|
||||
}
|
||||
}
|
||||
|
||||
p->totalParticles = totalParticles;
|
||||
p->pointer = totalParticles + 1;
|
||||
memcpy(p->particlePool, tempPool, totalParticles * sizeof(Particle));
|
||||
}
|
||||
|
@ -19,13 +19,12 @@ typedef struct {
|
||||
} Particle;
|
||||
|
||||
typedef struct {
|
||||
int numberOfParticles, totalParticles;
|
||||
int numParticlesInLine, removedParticles, totalParticles;
|
||||
double startTime, injectTimePeriod, writeTimePeriod;
|
||||
double lastInjectTime, lastUpdateTime, lastWriteTime;
|
||||
|
||||
int estimatedNumParticles;
|
||||
int numAllocatedParticles;
|
||||
|
||||
double dt;
|
||||
Particle* linSpaceLine;
|
||||
Particle* particlePool;
|
||||
|
||||
@ -35,12 +34,8 @@ typedef struct {
|
||||
} ParticleTracer;
|
||||
|
||||
extern void initParticleTracer(ParticleTracer*, Grid*, Parameter*);
|
||||
extern void injectParticles(ParticleTracer*);
|
||||
extern void advanceParticles(ParticleTracer*, double*, double*, int*, double);
|
||||
extern void freeParticles(ParticleTracer*);
|
||||
extern void writeParticles(ParticleTracer*);
|
||||
extern void printParticleTracerParameters(ParticleTracer*);
|
||||
extern void printParticles(ParticleTracer*);
|
||||
extern void trace(ParticleTracer*, double*, double*, int*, double);
|
||||
extern void compress(ParticleTracer*);
|
||||
extern void trace(ParticleTracer*, double*, double*, double, double);
|
||||
#endif
|
||||
|
@ -9,6 +9,6 @@
|
||||
|
||||
extern void initProgress(double);
|
||||
extern void printProgress(double);
|
||||
extern void stopProgress();
|
||||
extern void stopProgress(void);
|
||||
|
||||
#endif
|
||||
|
186
EnhancedSolver/2D-seq/src/solver-mg.c
Normal file
186
EnhancedSolver/2D-seq/src/solver-mg.c
Normal file
@ -0,0 +1,186 @@
|
||||
/*
|
||||
* Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
|
||||
* All rights reserved. This file is part of nusif-solver.
|
||||
* Use of this source code is governed by a MIT style
|
||||
* license that can be found in the LICENSE file.
|
||||
*/
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
|
||||
#include "allocate.h"
|
||||
#include "solver.h"
|
||||
#include "util.h"
|
||||
|
||||
#define FINEST_LEVEL 0
|
||||
#define COARSEST_LEVEL (s->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)
|
||||
{
|
||||
double* r = s->r[level + 1];
|
||||
double* old = s->r[level];
|
||||
|
||||
for (int j = 1; j < jmax + 1; j++) {
|
||||
for (int i = 1; i < imax + 1; i++) {
|
||||
R(i, j) = (OLD(2 * i - 1, 2 * j - 1) + OLD(2 * i, 2 * j - 1) * 2 +
|
||||
OLD(2 * i + 1, 2 * j - 1) + OLD(2 * i - 1, 2 * j) * 2 +
|
||||
OLD(2 * i, 2 * j) * 4 + OLD(2 * i + 1, 2 * j) * 2 +
|
||||
OLD(2 * i - 1, 2 * j + 1) + OLD(2 * i, 2 * j + 1) * 2 +
|
||||
OLD(2 * i + 1, 2 * j + 1)) /
|
||||
16.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
static void prolongate(Solver* s, int level, int imax, int jmax)
|
||||
{
|
||||
double* old = s->r[level + 1];
|
||||
double* e = s->r[level];
|
||||
|
||||
for (int j = 2; j < jmax + 1; j += 2) {
|
||||
for (int i = 2; i < imax + 1; i += 2) {
|
||||
E(i, j) = OLD(i / 2, j / 2);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
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);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
static void setBoundaryCondition(double* p, int imax, int jmax)
|
||||
{
|
||||
for (int i = 1; i < imax + 1; i++) {
|
||||
P(i, 0) = P(i, 1);
|
||||
P(i, jmax + 1) = P(i, jmax);
|
||||
}
|
||||
|
||||
for (int j = 1; j < jmax + 1; j++) {
|
||||
P(0, j) = P(1, j);
|
||||
P(imax + 1, j) = P(imax, j);
|
||||
}
|
||||
}
|
||||
|
||||
static double smooth(Solver* s, double* p, double* rhs, int level, int imax, int jmax)
|
||||
{
|
||||
double dx2 = s->grid->dx * s->grid->dx;
|
||||
double dy2 = s->grid->dy * s->grid->dy;
|
||||
double idx2 = 1.0 / dx2;
|
||||
double idy2 = 1.0 / dy2;
|
||||
double factor = s->omega * 0.5 * (dx2 * dy2) / (dx2 + dy2);
|
||||
double* r = s->r[level];
|
||||
double res = 1.0;
|
||||
int pass, jsw, isw;
|
||||
|
||||
jsw = 1;
|
||||
|
||||
for (pass = 0; pass < 2; pass++) {
|
||||
isw = jsw;
|
||||
|
||||
for (int j = 1; j < jmax + 1; j++) {
|
||||
for (int i = isw; i < imax + 1; i += 2) {
|
||||
|
||||
R(i, j) = RHS(i, j) -
|
||||
((P(i + 1, j) - 2.0 * P(i, j) + P(i - 1, j)) * idx2 +
|
||||
(P(i, j + 1) - 2.0 * P(i, j) + P(i, j - 1)) * idy2);
|
||||
|
||||
P(i, j) -= (factor * R(i, j));
|
||||
res += (R(i, j) * R(i, j));
|
||||
}
|
||||
isw = 3 - isw;
|
||||
}
|
||||
jsw = 3 - jsw;
|
||||
}
|
||||
|
||||
res = res / (double)(imax * jmax);
|
||||
return res;
|
||||
}
|
||||
|
||||
static double multiGrid(Solver* s, double* p, double* rhs, int level, int imax, int jmax)
|
||||
{
|
||||
double res = 0.0;
|
||||
|
||||
// coarsest level TODO: Use direct solver?
|
||||
if (level == COARSEST_LEVEL) {
|
||||
for (int i = 0; i < 5; i++) {
|
||||
smooth(s, p, rhs, level, imax, jmax);
|
||||
}
|
||||
return res;
|
||||
}
|
||||
|
||||
// pre-smoothing TODO: Make smoothing steps configurable?
|
||||
for (int i = 0; i < 5; i++) {
|
||||
smooth(s, p, rhs, level, imax, jmax);
|
||||
if (level == FINEST_LEVEL) setBoundaryCondition(p, imax, jmax);
|
||||
}
|
||||
|
||||
// restrict
|
||||
restrictMG(s, level, imax, jmax);
|
||||
|
||||
// MGSolver on residual and error.
|
||||
// TODO: What if there is a rest?
|
||||
multiGrid(s, s->e[level + 1], s->r[level + 1], level + 1, imax / 2, jmax / 2);
|
||||
|
||||
// prolongate
|
||||
prolongate(s, level, imax, jmax);
|
||||
|
||||
// correct p on finer level using residual
|
||||
correct(s, p, level, imax, jmax);
|
||||
if (level == FINEST_LEVEL) setBoundaryCondition(p, imax, jmax);
|
||||
|
||||
// post-smoothing
|
||||
for (int i = 0; i < 5; i++) {
|
||||
res = smooth(s, p, rhs, level, imax, jmax);
|
||||
if (level == FINEST_LEVEL) setBoundaryCondition(p, imax, jmax);
|
||||
}
|
||||
|
||||
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);
|
||||
|
||||
#ifdef VERBOSE
|
||||
printf("Residuum: %.6f\n", res);
|
||||
#endif
|
||||
}
|
128
EnhancedSolver/2D-seq/src/solver-sor.c
Normal file
128
EnhancedSolver/2D-seq/src/solver-sor.c
Normal file
@ -0,0 +1,128 @@
|
||||
/*
|
||||
* Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
|
||||
* All rights reserved. This file is part of nusif-solver.
|
||||
* Use of this source code is governed by a MIT style
|
||||
* license that can be found in the LICENSE file.
|
||||
*/
|
||||
#include "solver.h"
|
||||
#include "util.h"
|
||||
|
||||
void initSolver(Solver* s, Discretization* d, Parameter* p)
|
||||
{
|
||||
s->grid = &d->grid;
|
||||
s->itermax = p->itermax;
|
||||
s->eps = p->eps;
|
||||
s->omega = p->omg;
|
||||
}
|
||||
|
||||
void solveSOR(Solver* solver, double* p, double* rhs)
|
||||
{
|
||||
int imax = solver->grid->imax;
|
||||
int jmax = solver->grid->jmax;
|
||||
double eps = solver->eps;
|
||||
int itermax = solver->itermax;
|
||||
double dx2 = solver->grid->dx * solver->grid->dx;
|
||||
double dy2 = solver->grid->dy * solver->grid->dy;
|
||||
double idx2 = 1.0 / dx2;
|
||||
double idy2 = 1.0 / dy2;
|
||||
double factor = solver->omega * 0.5 * (dx2 * dy2) / (dx2 + dy2);
|
||||
double epssq = eps * eps;
|
||||
int it = 0;
|
||||
double res = 1.0;
|
||||
|
||||
while ((res >= epssq) && (it < itermax)) {
|
||||
res = 0.0;
|
||||
|
||||
for (int j = 1; j < jmax + 1; j++) {
|
||||
for (int i = 1; i < imax + 1; i++) {
|
||||
|
||||
double r = RHS(i, j) -
|
||||
((P(i + 1, j) - 2.0 * P(i, j) + P(i - 1, j)) * idx2 +
|
||||
(P(i, j + 1) - 2.0 * P(i, j) + P(i, j - 1)) * idy2);
|
||||
|
||||
P(i, j) -= (factor * r);
|
||||
res += (r * r);
|
||||
}
|
||||
}
|
||||
|
||||
for (int i = 1; i < imax + 1; i++) {
|
||||
P(i, 0) = P(i, 1);
|
||||
P(i, jmax + 1) = P(i, jmax);
|
||||
}
|
||||
|
||||
for (int j = 1; j < jmax + 1; j++) {
|
||||
P(0, j) = P(1, j);
|
||||
P(imax + 1, j) = P(imax, j);
|
||||
}
|
||||
|
||||
res = res / (double)(imax * jmax);
|
||||
#ifdef DEBUG
|
||||
printf("%d Residuum: %e\n", it, res);
|
||||
#endif
|
||||
it++;
|
||||
}
|
||||
|
||||
#ifdef VERBOSE
|
||||
printf("Solver took %d iterations to reach %f\n", it, sqrt(res));
|
||||
#endif
|
||||
}
|
||||
|
||||
void solve(Solver* solver, double* p, double* rhs)
|
||||
{
|
||||
int imax = solver->grid->imax;
|
||||
int jmax = solver->grid->jmax;
|
||||
double eps = solver->eps;
|
||||
int itermax = solver->itermax;
|
||||
double dx2 = solver->grid->dx * solver->grid->dx;
|
||||
double dy2 = solver->grid->dy * solver->grid->dy;
|
||||
double idx2 = 1.0 / dx2;
|
||||
double idy2 = 1.0 / dy2;
|
||||
double factor = solver->omega * 0.5 * (dx2 * dy2) / (dx2 + dy2);
|
||||
double epssq = eps * eps;
|
||||
int it = 0;
|
||||
double res = 1.0;
|
||||
int pass, jsw, isw;
|
||||
|
||||
while ((res >= epssq) && (it < itermax)) {
|
||||
res = 0.0;
|
||||
jsw = 1;
|
||||
|
||||
for (pass = 0; pass < 2; pass++) {
|
||||
isw = jsw;
|
||||
|
||||
for (int j = 1; j < jmax + 1; j++) {
|
||||
for (int i = isw; i < imax + 1; i += 2) {
|
||||
|
||||
double r = RHS(i, j) -
|
||||
((P(i + 1, j) - 2.0 * P(i, j) + P(i - 1, j)) * idx2 +
|
||||
(P(i, j + 1) - 2.0 * P(i, j) + P(i, j - 1)) * idy2);
|
||||
|
||||
P(i, j) -= (factor * r);
|
||||
res += (r * r);
|
||||
}
|
||||
isw = 3 - isw;
|
||||
}
|
||||
jsw = 3 - jsw;
|
||||
}
|
||||
|
||||
for (int i = 1; i < imax + 1; i++) {
|
||||
P(i, 0) = P(i, 1);
|
||||
P(i, jmax + 1) = P(i, jmax);
|
||||
}
|
||||
|
||||
for (int j = 1; j < jmax + 1; j++) {
|
||||
P(0, j) = P(1, j);
|
||||
P(imax + 1, j) = P(imax, j);
|
||||
}
|
||||
|
||||
res = res / (double)(imax * jmax);
|
||||
#ifdef DEBUG
|
||||
printf("%d Residuum: %e\n", it, res);
|
||||
#endif
|
||||
it++;
|
||||
}
|
||||
|
||||
#ifdef VERBOSE
|
||||
printf("Solver took %d iterations to reach %f\n", it, sqrt(res));
|
||||
#endif
|
||||
}
|
@ -1,900 +0,0 @@
|
||||
/*
|
||||
* Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
|
||||
* All rights reserved. This file is part of nusif-solver.
|
||||
* Use of this source code is governed by a MIT style
|
||||
* license that can be found in the LICENSE file.
|
||||
*/
|
||||
#include <float.h>
|
||||
#include <math.h>
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
|
||||
#include "allocate.h"
|
||||
#include "parameter.h"
|
||||
#include "solver.h"
|
||||
#include "util.h"
|
||||
|
||||
#define P(i, j) p[(j) * (imax + 2) + (i)]
|
||||
#define F(i, j) f[(j) * (imax + 2) + (i)]
|
||||
#define G(i, j) g[(j) * (imax + 2) + (i)]
|
||||
#define U(i, j) u[(j) * (imax + 2) + (i)]
|
||||
#define V(i, j) v[(j) * (imax + 2) + (i)]
|
||||
#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 oldR(i, j) oldr[(j) * (imax + 2) + (i)]
|
||||
#define oldE(i, j) olde[(j) * (imax + 2) + (i)]
|
||||
#define RHS(i, j) rhs[(j) * (imax + 2) + (i)]
|
||||
|
||||
static double distance(double i, double j, double iCenter, double jCenter)
|
||||
{
|
||||
return sqrt(pow(iCenter - i, 2) + pow(jCenter - j, 2) * 1.0);
|
||||
}
|
||||
|
||||
void print(Solver* solver, double* grid)
|
||||
{
|
||||
int imax = solver->grid.imax;
|
||||
int jmax = solver->grid.jmax;
|
||||
|
||||
for (int j = 0; j < jmax + 2; j++) {
|
||||
printf("%02d: ", j);
|
||||
for (int i = 0; i < imax + 2; i++) {
|
||||
printf("%3.2f ", grid[j * (imax + 2) + i]);
|
||||
}
|
||||
printf("\n");
|
||||
}
|
||||
fflush(stdout);
|
||||
}
|
||||
|
||||
void printGrid(Solver* solver, int* grid)
|
||||
{
|
||||
int imax = solver->grid.imax;
|
||||
int jmax = solver->grid.jmax;
|
||||
|
||||
for (int j = 0; j < jmax + 2; j++) {
|
||||
printf("%02d: ", j);
|
||||
for (int i = 0; i < imax + 2; i++) {
|
||||
printf("%2d ", grid[j * (imax + 2) + i]);
|
||||
}
|
||||
printf("\n");
|
||||
}
|
||||
fflush(stdout);
|
||||
}
|
||||
|
||||
static void printConfig(Solver* solver)
|
||||
{
|
||||
printf("Parameters for #%s#\n", solver->problem);
|
||||
printf("Boundary conditions Left:%d Right:%d Bottom:%d Top:%d\n",
|
||||
solver->bcLeft,
|
||||
solver->bcRight,
|
||||
solver->bcBottom,
|
||||
solver->bcTop);
|
||||
printf("\tReynolds number: %.2f\n", solver->re);
|
||||
printf("\tGx Gy: %.2f %.2f\n", solver->gx, solver->gy);
|
||||
printf("Geometry data:\n");
|
||||
printf("\tDomain box size (x, y): %.2f, %.2f\n",
|
||||
solver->grid.xlength,
|
||||
solver->grid.ylength);
|
||||
printf("\tCells (x, y): %d, %d\n", solver->grid.imax, solver->grid.jmax);
|
||||
printf("Timestep parameters:\n");
|
||||
printf("\tDefault stepsize: %.2f, Final time %.2f\n", solver->dt, solver->te);
|
||||
printf("\tdt bound: %.6f\n", solver->dtBound);
|
||||
printf("\tTau factor: %.2f\n", solver->tau);
|
||||
printf("Iterative solver parameters:\n");
|
||||
printf("\tMax iterations: %d\n", solver->itermax);
|
||||
printf("\tepsilon (stopping tolerance) : %f\n", solver->eps);
|
||||
printf("\tgamma factor: %f\n", solver->gamma);
|
||||
printf("\tomega (SOR relaxation): %f\n", solver->omega);
|
||||
}
|
||||
|
||||
void initSolver(Solver* solver, Parameter* params)
|
||||
{
|
||||
solver->problem = params->name;
|
||||
solver->bcLeft = params->bcLeft;
|
||||
solver->bcRight = params->bcRight;
|
||||
solver->bcBottom = params->bcBottom;
|
||||
solver->bcTop = params->bcTop;
|
||||
solver->grid.imax = params->imax;
|
||||
solver->grid.jmax = params->jmax;
|
||||
solver->grid.xlength = params->xlength;
|
||||
solver->grid.ylength = params->ylength;
|
||||
solver->grid.dx = params->xlength / params->imax;
|
||||
solver->grid.dy = params->ylength / params->jmax;
|
||||
solver->eps = params->eps;
|
||||
solver->omega = params->omg;
|
||||
solver->itermax = params->itermax;
|
||||
solver->re = params->re;
|
||||
solver->gx = params->gx;
|
||||
solver->gy = params->gy;
|
||||
solver->dt = params->dt;
|
||||
solver->te = params->te;
|
||||
solver->tau = params->tau;
|
||||
solver->gamma = params->gamma;
|
||||
solver->rho = params->rho;
|
||||
solver->levels = params->levels;
|
||||
solver->currentlevel = 0;
|
||||
|
||||
int imax = solver->grid.imax;
|
||||
int jmax = solver->grid.jmax;
|
||||
int levels = solver->levels;
|
||||
|
||||
size_t size_level = levels * (imax + 2) * (jmax + 2) * sizeof(double);
|
||||
|
||||
size_t size = (imax + 2) * (jmax + 2) * sizeof(double);
|
||||
solver->u = allocate(64, size);
|
||||
solver->v = allocate(64, size);
|
||||
solver->s = allocate(64, size);
|
||||
solver->p = allocate(64, size);
|
||||
solver->rhs = allocate(64, size);
|
||||
solver->f = allocate(64, size);
|
||||
solver->g = allocate(64, size);
|
||||
|
||||
solver->r = malloc(levels * sizeof(double*));
|
||||
solver->e = malloc(levels * sizeof(double*));
|
||||
|
||||
for (int j = 0; j < levels; ++j) {
|
||||
|
||||
solver->r[j] = allocate(64, size);
|
||||
solver->e[j] = allocate(64, size);
|
||||
}
|
||||
|
||||
for (int i = 0; i < (imax + 2) * (jmax + 2); i++) {
|
||||
|
||||
solver->u[i] = params->u_init;
|
||||
solver->v[i] = params->v_init;
|
||||
solver->p[i] = params->p_init;
|
||||
solver->rhs[i] = 0.0;
|
||||
solver->f[i] = 0.0;
|
||||
solver->g[i] = 0.0;
|
||||
solver->s[i] = NONE;
|
||||
for (int j = 0; j < levels; ++j) {
|
||||
|
||||
solver->r[j][i] = 0.0;
|
||||
solver->e[j][i] = 0.0;
|
||||
}
|
||||
}
|
||||
|
||||
double dx = solver->grid.dx;
|
||||
double dy = solver->grid.dy;
|
||||
double invSqrSum = 1.0 / (dx * dx) + 1.0 / (dy * dy);
|
||||
solver->dtBound = 0.5 * solver->re * 1.0 / invSqrSum;
|
||||
|
||||
double xCenter = 0, yCenter = 0, radius = 0;
|
||||
double x1 = 0, x2 = 0, y1 = 0, y2 = 0;
|
||||
|
||||
int* s = solver->s;
|
||||
|
||||
switch (params->shape) {
|
||||
case NOSHAPE:
|
||||
break;
|
||||
case RECT:
|
||||
x1 = params->xCenter - params->xRectLength / 2;
|
||||
x2 = params->xCenter + params->xRectLength / 2;
|
||||
y1 = params->yCenter - params->yRectLength / 2;
|
||||
y2 = params->yCenter + params->yRectLength / 2;
|
||||
|
||||
for (int j = 1; j < jmax + 1; j++) {
|
||||
for (int i = 1; i < imax + 1; i++) {
|
||||
if ((x1 <= (i * dx)) && ((i * dx) <= x2) && (y1 <= (j * dy)) &&
|
||||
((j * dy) <= y2)) {
|
||||
S(i, j) = LOCAL;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
break;
|
||||
case CIRCLE:
|
||||
xCenter = params->xCenter;
|
||||
yCenter = params->yCenter;
|
||||
radius = params->circleRadius;
|
||||
|
||||
for (int j = 1; j < jmax + 1; j++) {
|
||||
for (int i = 1; i < imax + 1; i++) {
|
||||
if (distance((i * dx), (j * dy), xCenter, yCenter) <= radius) {
|
||||
S(i, j) = LOCAL;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
break;
|
||||
default:
|
||||
break;
|
||||
}
|
||||
|
||||
if (params->shape != NOSHAPE) {
|
||||
for (int j = 1; j < jmax + 1; j++) {
|
||||
for (int i = 1; i < imax + 1; i++) {
|
||||
|
||||
if (S(i, j - 1) == NONE && S(i, j + 1) == LOCAL && S(i, j) == LOCAL)
|
||||
S(i, j) = BOTTOM; // TOP
|
||||
if (S(i - 1, j) == NONE && S(i + 1, j) == LOCAL && S(i, j) == LOCAL)
|
||||
S(i, j) = LEFT; // LEFT
|
||||
if (S(i + 1, j) == NONE && S(i - 1, j) == LOCAL && S(i, j) == LOCAL)
|
||||
S(i, j) = RIGHT; // RIGHT
|
||||
if (S(i, j + 1) == NONE && S(i, j - 1) == LOCAL && S(i, j) == LOCAL)
|
||||
S(i, j) = TOP; // BOTTOM
|
||||
if (S(i - 1, j - 1) == NONE && S(i, j - 1) == NONE &&
|
||||
S(i - 1, j) == NONE && S(i + 1, j + 1) == LOCAL &&
|
||||
(S(i, j) == LOCAL || S(i, j) == LEFT || S(i, j) == BOTTOM))
|
||||
S(i, j) = BOTTOMLEFT; // TOPLEFT
|
||||
if (S(i + 1, j - 1) == NONE && S(i, j - 1) == NONE &&
|
||||
S(i + 1, j) == NONE && S(i - 1, j + 1) == LOCAL &&
|
||||
(S(i, j) == LOCAL || S(i, j) == RIGHT || S(i, j) == BOTTOM))
|
||||
S(i, j) = BOTTOMRIGHT; // TOPRIGHT
|
||||
if (S(i - 1, j + 1) == NONE && S(i - 1, j) == NONE &&
|
||||
S(i, j + 1) == NONE && S(i + 1, j - 1) == LOCAL &&
|
||||
(S(i, j) == LOCAL || S(i, j) == LEFT || S(i, j) == TOP))
|
||||
S(i, j) = TOPLEFT; // BOTTOMLEFT
|
||||
if (S(i + 1, j + 1) == NONE && S(i + 1, j) == NONE &&
|
||||
S(i, j + 1) == NONE && S(i - 1, j - 1) == LOCAL &&
|
||||
(S(i, j) == LOCAL || S(i, j) == RIGHT || S(i, j) == TOP))
|
||||
S(i, j) = TOPRIGHT; // BOTTOMRIGHT
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#ifdef VERBOSE
|
||||
printConfig(solver);
|
||||
#endif
|
||||
}
|
||||
|
||||
static double maxElement(Solver* solver, double* m)
|
||||
{
|
||||
int size = (solver->grid.imax + 2) * (solver->grid.jmax + 2);
|
||||
double maxval = DBL_MIN;
|
||||
|
||||
for (int i = 0; i < size; i++) {
|
||||
maxval = MAX(maxval, fabs(m[i]));
|
||||
}
|
||||
|
||||
return maxval;
|
||||
}
|
||||
|
||||
void computeRHS(Solver* solver)
|
||||
{
|
||||
int imax = solver->grid.imax;
|
||||
int jmax = solver->grid.jmax;
|
||||
double idx = 1.0 / solver->grid.dx;
|
||||
double idy = 1.0 / solver->grid.dy;
|
||||
double idt = 1.0 / solver->dt;
|
||||
double* rhs = solver->rhs;
|
||||
double* f = solver->f;
|
||||
double* g = solver->g;
|
||||
int* s = solver->s;
|
||||
|
||||
for (int j = 1; j < jmax + 1; j++) {
|
||||
for (int i = 1; i < imax + 1; i++) {
|
||||
RHS(i, j) = idt *
|
||||
((F(i, j) - F(i - 1, j)) * idx + (G(i, j) - G(i, j - 1)) * idy);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void normalizePressure(Solver* solver)
|
||||
{
|
||||
int size = (solver->grid.imax + 2) * (solver->grid.jmax + 2);
|
||||
double* p = solver->p;
|
||||
double avgP = 0.0;
|
||||
|
||||
for (int i = 0; i < size; i++) {
|
||||
avgP += p[i];
|
||||
}
|
||||
avgP /= size;
|
||||
|
||||
for (int i = 0; i < size; i++) {
|
||||
p[i] = p[i] - avgP;
|
||||
}
|
||||
}
|
||||
|
||||
void computeTimestep(Solver* solver)
|
||||
{
|
||||
double dt = solver->dtBound;
|
||||
double dx = solver->grid.dx;
|
||||
double dy = solver->grid.dy;
|
||||
double umax = maxElement(solver, solver->u);
|
||||
double vmax = maxElement(solver, solver->v);
|
||||
|
||||
if (umax > 0) {
|
||||
dt = (dt > dx / umax) ? dx / umax : dt;
|
||||
}
|
||||
if (vmax > 0) {
|
||||
dt = (dt > dy / vmax) ? dy / vmax : dt;
|
||||
}
|
||||
|
||||
solver->dt = dt * solver->tau;
|
||||
}
|
||||
|
||||
void setBoundaryConditions(Solver* solver)
|
||||
{
|
||||
int imax = solver->grid.imax;
|
||||
int jmax = solver->grid.jmax;
|
||||
double* u = solver->u;
|
||||
double* v = solver->v;
|
||||
|
||||
// Left boundary
|
||||
switch (solver->bcLeft) {
|
||||
case NOSLIP:
|
||||
for (int j = 1; j < jmax + 1; j++) {
|
||||
U(0, j) = 0.0;
|
||||
V(0, j) = -V(1, j);
|
||||
}
|
||||
break;
|
||||
case SLIP:
|
||||
for (int j = 1; j < jmax + 1; j++) {
|
||||
U(0, j) = 0.0;
|
||||
V(0, j) = V(1, j);
|
||||
}
|
||||
break;
|
||||
case OUTFLOW:
|
||||
for (int j = 1; j < jmax + 1; j++) {
|
||||
U(0, j) = U(1, j);
|
||||
V(0, j) = V(1, j);
|
||||
}
|
||||
break;
|
||||
case PERIODIC:
|
||||
break;
|
||||
}
|
||||
|
||||
// Right boundary
|
||||
switch (solver->bcRight) {
|
||||
case NOSLIP:
|
||||
for (int j = 1; j < jmax + 1; j++) {
|
||||
U(imax, j) = 0.0;
|
||||
V(imax + 1, j) = -V(imax, j);
|
||||
}
|
||||
break;
|
||||
case SLIP:
|
||||
for (int j = 1; j < jmax + 1; j++) {
|
||||
U(imax, j) = 0.0;
|
||||
V(imax + 1, j) = V(imax, j);
|
||||
}
|
||||
break;
|
||||
case OUTFLOW:
|
||||
for (int j = 1; j < jmax + 1; j++) {
|
||||
U(imax, j) = U(imax - 1, j);
|
||||
V(imax + 1, j) = V(imax, j);
|
||||
}
|
||||
break;
|
||||
case PERIODIC:
|
||||
break;
|
||||
}
|
||||
|
||||
// Bottom boundary
|
||||
switch (solver->bcBottom) {
|
||||
case NOSLIP:
|
||||
for (int i = 1; i < imax + 1; i++) {
|
||||
V(i, 0) = 0.0;
|
||||
U(i, 0) = -U(i, 1);
|
||||
}
|
||||
break;
|
||||
case SLIP:
|
||||
for (int i = 1; i < imax + 1; i++) {
|
||||
V(i, 0) = 0.0;
|
||||
U(i, 0) = U(i, 1);
|
||||
}
|
||||
break;
|
||||
case OUTFLOW:
|
||||
for (int i = 1; i < imax + 1; i++) {
|
||||
U(i, 0) = U(i, 1);
|
||||
V(i, 0) = V(i, 1);
|
||||
}
|
||||
break;
|
||||
case PERIODIC:
|
||||
break;
|
||||
}
|
||||
|
||||
// Top boundary
|
||||
switch (solver->bcTop) {
|
||||
case NOSLIP:
|
||||
for (int i = 1; i < imax + 1; i++) {
|
||||
V(i, jmax) = 0.0;
|
||||
U(i, jmax + 1) = -U(i, jmax);
|
||||
}
|
||||
break;
|
||||
case SLIP:
|
||||
for (int i = 1; i < imax + 1; i++) {
|
||||
V(i, jmax) = 0.0;
|
||||
U(i, jmax + 1) = U(i, jmax);
|
||||
}
|
||||
break;
|
||||
case OUTFLOW:
|
||||
for (int i = 1; i < imax + 1; i++) {
|
||||
U(i, jmax + 1) = U(i, jmax);
|
||||
V(i, jmax) = V(i, jmax - 1);
|
||||
}
|
||||
break;
|
||||
case PERIODIC:
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
void setSpecialBoundaryCondition(Solver* solver)
|
||||
{
|
||||
int imax = solver->grid.imax;
|
||||
int jmax = solver->grid.jmax;
|
||||
double mDy = solver->grid.dy;
|
||||
double* u = solver->u;
|
||||
int* s = solver->s;
|
||||
|
||||
if (strcmp(solver->problem, "dcavity") == 0) {
|
||||
for (int i = 1; i < imax; i++) {
|
||||
U(i, jmax + 1) = 2.0 - U(i, jmax);
|
||||
}
|
||||
} else if (strcmp(solver->problem, "canal") == 0) {
|
||||
double ylength = solver->grid.ylength;
|
||||
double y;
|
||||
|
||||
for (int j = 1; j < jmax + 1; j++) {
|
||||
y = mDy * (j - 0.5);
|
||||
U(0, j) = y * (ylength - y) * 4.0 / (ylength * ylength);
|
||||
}
|
||||
} else if (strcmp(solver->problem, "backstep") == 0) {
|
||||
for (int j = 1; j < jmax + 1; j++) {
|
||||
if (S(0, j) == NONE) U(0, j) = 1.0;
|
||||
}
|
||||
} else if (strcmp(solver->problem, "karman") == 0) {
|
||||
for (int j = 1; j < jmax + 1; j++) {
|
||||
U(0, j) = 1.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void setObjectBoundaryCondition(Solver* solver)
|
||||
{
|
||||
int imax = solver->grid.imax;
|
||||
int jmax = solver->grid.jmax;
|
||||
double* u = solver->u;
|
||||
double* v = solver->v;
|
||||
int* s = solver->s;
|
||||
|
||||
for (int j = 1; j < jmax + 1; j++) {
|
||||
for (int i = 1; i < imax + 1; i++) {
|
||||
switch (S(i, j)) {
|
||||
case TOP:
|
||||
U(i, j) = -U(i, j + 1);
|
||||
U(i - 1, j) = -U(i - 1, j + 1);
|
||||
V(i, j) = 0.0;
|
||||
break;
|
||||
case BOTTOM:
|
||||
U(i, j) = -U(i, j - 1);
|
||||
U(i - 1, j) = -U(i - 1, j - 1);
|
||||
V(i, j) = 0.0;
|
||||
break;
|
||||
case LEFT:
|
||||
U(i - 1, j) = 0.0;
|
||||
V(i, j) = -V(i - 1, j);
|
||||
V(i, j - 1) = -V(i - 1, j - 1);
|
||||
break;
|
||||
case RIGHT:
|
||||
U(i, j) = 0.0;
|
||||
V(i, j) = -V(i + 1, j);
|
||||
V(i, j - 1) = -V(i + 1, j - 1);
|
||||
break;
|
||||
case TOPLEFT:
|
||||
U(i, j) = -U(i, j + 1);
|
||||
U(i - 1, j) = 0.0;
|
||||
V(i, j) = 0.0;
|
||||
V(i, j - 1) = -V(i - 1, j - 1);
|
||||
break;
|
||||
case TOPRIGHT:
|
||||
U(i, j) = 0.0;
|
||||
U(i - 1, j) = -U(i - 1, j + 1);
|
||||
V(i, j) = 0.0;
|
||||
V(i, j - 1) = -V(i + 1, j - 1);
|
||||
break;
|
||||
case BOTTOMLEFT:
|
||||
U(i, j) = -U(i, j - 1);
|
||||
U(i - 1, j) = 0.0;
|
||||
V(i, j) = -V(i - 1, j);
|
||||
V(i, j - 1) = 0.0;
|
||||
break;
|
||||
case BOTTOMRIGHT:
|
||||
U(i, j) = 0.0;
|
||||
U(i - 1, j) = -U(i - 1, j - 1);
|
||||
V(i, j) = -V(i, j + 1);
|
||||
V(i, j - 1) = 0.0;
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void computeFG(Solver* solver)
|
||||
{
|
||||
double* u = solver->u;
|
||||
double* v = solver->v;
|
||||
double* f = solver->f;
|
||||
double* g = solver->g;
|
||||
int* s = solver->s;
|
||||
int imax = solver->grid.imax;
|
||||
int jmax = solver->grid.jmax;
|
||||
double gx = solver->gx;
|
||||
double gy = solver->gy;
|
||||
double gamma = solver->gamma;
|
||||
double dt = solver->dt;
|
||||
double inverseRe = 1.0 / solver->re;
|
||||
double inverseDx = 1.0 / solver->grid.dx;
|
||||
double inverseDy = 1.0 / solver->grid.dy;
|
||||
double du2dx, dv2dy, duvdx, duvdy;
|
||||
double du2dx2, du2dy2, dv2dx2, dv2dy2;
|
||||
|
||||
for (int j = 1; j < jmax + 1; j++) {
|
||||
for (int i = 1; i < imax + 1; i++) {
|
||||
if (S(i, j) == NONE) {
|
||||
du2dx = inverseDx * 0.25 *
|
||||
((U(i, j) + U(i + 1, j)) * (U(i, j) + U(i + 1, j)) -
|
||||
(U(i, j) + U(i - 1, j)) * (U(i, j) + U(i - 1, j))) +
|
||||
gamma * inverseDx * 0.25 *
|
||||
(fabs(U(i, j) + U(i + 1, j)) * (U(i, j) - U(i + 1, j)) +
|
||||
fabs(U(i, j) + U(i - 1, j)) * (U(i, j) - U(i - 1, j)));
|
||||
|
||||
duvdy = inverseDy * 0.25 *
|
||||
((V(i, j) + V(i + 1, j)) * (U(i, j) + U(i, j + 1)) -
|
||||
(V(i, j - 1) + V(i + 1, j - 1)) *
|
||||
(U(i, j) + U(i, j - 1))) +
|
||||
gamma * inverseDy * 0.25 *
|
||||
(fabs(V(i, j) + V(i + 1, j)) * (U(i, j) - U(i, j + 1)) +
|
||||
fabs(V(i, j - 1) + V(i + 1, j - 1)) *
|
||||
(U(i, j) - U(i, j - 1)));
|
||||
|
||||
du2dx2 = inverseDx * inverseDx *
|
||||
(U(i + 1, j) - 2.0 * U(i, j) + U(i - 1, j));
|
||||
du2dy2 = inverseDy * inverseDy *
|
||||
(U(i, j + 1) - 2.0 * U(i, j) + U(i, j - 1));
|
||||
F(i, j) = U(i, j) +
|
||||
dt * (inverseRe * (du2dx2 + du2dy2) - du2dx - duvdy + gx);
|
||||
|
||||
duvdx = inverseDx * 0.25 *
|
||||
((U(i, j) + U(i, j + 1)) * (V(i, j) + V(i + 1, j)) -
|
||||
(U(i - 1, j) + U(i - 1, j + 1)) *
|
||||
(V(i, j) + V(i - 1, j))) +
|
||||
gamma * inverseDx * 0.25 *
|
||||
(fabs(U(i, j) + U(i, j + 1)) * (V(i, j) - V(i + 1, j)) +
|
||||
fabs(U(i - 1, j) + U(i - 1, j + 1)) *
|
||||
(V(i, j) - V(i - 1, j)));
|
||||
|
||||
dv2dy = inverseDy * 0.25 *
|
||||
((V(i, j) + V(i, j + 1)) * (V(i, j) + V(i, j + 1)) -
|
||||
(V(i, j) + V(i, j - 1)) * (V(i, j) + V(i, j - 1))) +
|
||||
gamma * inverseDy * 0.25 *
|
||||
(fabs(V(i, j) + V(i, j + 1)) * (V(i, j) - V(i, j + 1)) +
|
||||
fabs(V(i, j) + V(i, j - 1)) * (V(i, j) - V(i, j - 1)));
|
||||
|
||||
dv2dx2 = inverseDx * inverseDx *
|
||||
(V(i + 1, j) - 2.0 * V(i, j) + V(i - 1, j));
|
||||
dv2dy2 = inverseDy * inverseDy *
|
||||
(V(i, j + 1) - 2.0 * V(i, j) + V(i, j - 1));
|
||||
G(i, j) = V(i, j) +
|
||||
dt * (inverseRe * (dv2dx2 + dv2dy2) - duvdx - dv2dy + gy);
|
||||
} else {
|
||||
switch (S(i, j)) {
|
||||
case TOP:
|
||||
G(i, j) = V(i, j);
|
||||
break;
|
||||
case BOTTOM:
|
||||
G(i, j - 1) = V(i, j - 1);
|
||||
break;
|
||||
case LEFT:
|
||||
F(i - 1, j) = U(i - 1, j);
|
||||
break;
|
||||
case RIGHT:
|
||||
F(i, j) = U(i, j);
|
||||
break;
|
||||
case TOPLEFT:
|
||||
F(i - 1, j) = U(i - 1, j);
|
||||
G(i, j) = V(i, j);
|
||||
break;
|
||||
case TOPRIGHT:
|
||||
F(i, j) = U(i, j);
|
||||
G(i, j) = V(i, j);
|
||||
break;
|
||||
case BOTTOMLEFT:
|
||||
F(i - 1, j) = U(i - 1, j);
|
||||
G(i, j - 1) = V(i, j - 1);
|
||||
break;
|
||||
case BOTTOMRIGHT:
|
||||
F(i, j) = U(i, j);
|
||||
G(i, j - 1) = V(i, j - 1);
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------- boundary of F --------------------------- */
|
||||
for (int j = 1; j < jmax + 1; j++) {
|
||||
F(0, j) = U(0, j);
|
||||
F(imax, j) = U(imax, j);
|
||||
}
|
||||
|
||||
/* ---------------------- boundary of G --------------------------- */
|
||||
for (int i = 1; i < imax + 1; i++) {
|
||||
G(i, 0) = V(i, 0);
|
||||
G(i, jmax) = V(i, jmax);
|
||||
}
|
||||
}
|
||||
|
||||
void adaptUV(Solver* solver)
|
||||
{
|
||||
int imax = solver->grid.imax;
|
||||
int jmax = solver->grid.jmax;
|
||||
double* p = solver->p;
|
||||
double* u = solver->u;
|
||||
double* v = solver->v;
|
||||
int* s = solver->s;
|
||||
double* f = solver->f;
|
||||
double* g = solver->g;
|
||||
double factorX = solver->dt / solver->grid.dx;
|
||||
double factorY = solver->dt / solver->grid.dy;
|
||||
|
||||
double val = 0;
|
||||
|
||||
for (int j = 1; j < jmax + 1; j++) {
|
||||
for (int i = 1; i < imax + 1; i++) {
|
||||
U(i, j) = F(i, j) - (P(i + 1, j) - P(i, j)) * factorX;
|
||||
V(i, j) = G(i, j) - (P(i, j + 1) - P(i, j)) * factorY;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
double smoothRB(Solver* solver)
|
||||
{
|
||||
int imax = solver->grid.imax;
|
||||
int jmax = solver->grid.jmax;
|
||||
double eps = solver->eps;
|
||||
int itermax = solver->itermax;
|
||||
double dx2 = solver->grid.dx * solver->grid.dx;
|
||||
double dy2 = solver->grid.dy * solver->grid.dy;
|
||||
double idx2 = 1.0 / dx2;
|
||||
double idy2 = 1.0 / dy2;
|
||||
double factor = solver->omega * 0.5 * (dx2 * dy2) / (dx2 + dy2);
|
||||
double* p = solver->p;
|
||||
double* r = solver->r[solver->currentlevel];
|
||||
double* rhs = solver->rhs;
|
||||
double epssq = eps * eps;
|
||||
int it = 0;
|
||||
double res = 1.0;
|
||||
int pass, jsw, isw;
|
||||
|
||||
jsw = 1;
|
||||
|
||||
for (pass = 0; pass < 2; pass++) {
|
||||
isw = jsw;
|
||||
|
||||
for (int j = 1; j < jmax + 1; j++) {
|
||||
for (int i = isw; i < imax + 1; i += 2) {
|
||||
|
||||
R(i, j) = RHS(i, j) -
|
||||
((P(i + 1, j) - 2.0 * P(i, j) + P(i - 1, j)) * idx2 +
|
||||
(P(i, j + 1) - 2.0 * P(i, j) + P(i, j - 1)) * idy2);
|
||||
|
||||
P(i, j) -= (factor * R(i, j));
|
||||
res += (R(i, j) * R(i, j));
|
||||
}
|
||||
isw = 3 - isw;
|
||||
}
|
||||
jsw = 3 - jsw;
|
||||
}
|
||||
|
||||
res = res / (double)(imax * jmax);
|
||||
return res;
|
||||
}
|
||||
|
||||
void multiGrid(Solver* solver)
|
||||
{
|
||||
double res = 0.0;
|
||||
int imax = solver->grid.imax;
|
||||
int jmax = solver->grid.jmax;
|
||||
if (solver->currentlevel == (solver->levels - 1)) {
|
||||
for (int i = 0; i < 5; i++) {
|
||||
smoothRB(solver);
|
||||
}
|
||||
return;
|
||||
}
|
||||
|
||||
for (int i = 0; i < 5; i++) {
|
||||
smoothRB(solver);
|
||||
if (solver->currentlevel == 0) {
|
||||
|
||||
double* p = solver->p;
|
||||
|
||||
for (int i = 1; i < imax + 1; i++) {
|
||||
P(i, 0) = P(i, 1);
|
||||
P(i, jmax + 1) = P(i, jmax);
|
||||
}
|
||||
|
||||
for (int j = 1; j < jmax + 1; j++) {
|
||||
P(0, j) = P(1, j);
|
||||
P(imax + 1, j) = P(imax, j);
|
||||
}
|
||||
}
|
||||
}
|
||||
Solver coarseSolver = copySolver(solver);
|
||||
|
||||
// restrict
|
||||
restrictMG(solver);
|
||||
|
||||
coarseSolver.p = solver->e[coarseSolver.currentlevel];
|
||||
coarseSolver.rhs = solver->r[coarseSolver.currentlevel];
|
||||
coarseSolver.grid.imax /= 2;
|
||||
coarseSolver.grid.jmax /= 2;
|
||||
|
||||
// MGSolver on residual and error.
|
||||
multiGrid(&coarseSolver);
|
||||
|
||||
// prolongate
|
||||
prolongate(solver);
|
||||
|
||||
// correct p on finest level using residual
|
||||
correct(solver);
|
||||
|
||||
if (solver->currentlevel == 0) {
|
||||
|
||||
double* p = solver->p;
|
||||
|
||||
for (int i = 1; i < imax + 1; i++) {
|
||||
P(i, 0) = P(i, 1);
|
||||
P(i, jmax + 1) = P(i, jmax);
|
||||
}
|
||||
|
||||
for (int j = 1; j < jmax + 1; j++) {
|
||||
P(0, j) = P(1, j);
|
||||
P(imax + 1, j) = P(imax, j);
|
||||
}
|
||||
}
|
||||
|
||||
for (int i = 0; i < 5; i++) {
|
||||
res = smoothRB(solver);
|
||||
if (solver->currentlevel == 0) {
|
||||
|
||||
double* p = solver->p;
|
||||
|
||||
for (int i = 1; i < imax + 1; i++) {
|
||||
P(i, 0) = P(i, 1);
|
||||
P(i, jmax + 1) = P(i, jmax);
|
||||
}
|
||||
|
||||
for (int j = 1; j < jmax + 1; j++) {
|
||||
P(0, j) = P(1, j);
|
||||
P(imax + 1, j) = P(imax, j);
|
||||
}
|
||||
}
|
||||
}
|
||||
#ifdef VERBOSE
|
||||
if (solver->currentlevel == 0) {
|
||||
printf("Residuum: %.6f\n", res);
|
||||
}
|
||||
#endif
|
||||
}
|
||||
|
||||
void restrictMG(Solver* solver)
|
||||
{
|
||||
int imax = solver->grid.imax;
|
||||
int jmax = solver->grid.jmax;
|
||||
double* r = solver->r[solver->currentlevel + 1];
|
||||
double* oldr = solver->r[solver->currentlevel];
|
||||
|
||||
for (int j = 1; j < jmax + 1; j++) {
|
||||
for (int i = 1; i < imax + 1; ++i) {
|
||||
R(i, j) = (oldR(2 * i - 1, 2 * j - 1) + oldR(2 * i, 2 * j - 1) * 2 +
|
||||
oldR(2 * i + 1, 2 * j - 1) + oldR(2 * i - 1, 2 * j) * 2 +
|
||||
oldR(2 * i, 2 * j) * 4 + oldR(2 * i + 1, 2 * j) * 2 +
|
||||
oldR(2 * i - 1, 2 * j + 1) + oldR(2 * i, 2 * j + 1) * 2 +
|
||||
oldR(2 * i + 1, 2 * j + 1)) /
|
||||
16.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void prolongate(Solver* solver)
|
||||
{
|
||||
int imax = solver->grid.imax;
|
||||
int jmax = solver->grid.jmax;
|
||||
double* olde = solver->r[solver->currentlevel + 1];
|
||||
double* e = solver->r[solver->currentlevel];
|
||||
|
||||
for (int j = 2; j < jmax + 1; j += 2) {
|
||||
for (int i = 2; i < imax + 1; i += 2) {
|
||||
E(i, j) = oldE(i / 2, j / 2);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void correct(Solver* solver)
|
||||
{
|
||||
int imax = solver->grid.imax;
|
||||
int jmax = solver->grid.jmax;
|
||||
double* p = solver->p;
|
||||
double* e = solver->e[solver->currentlevel];
|
||||
for (int j = 1; j < jmax + 1; ++j) {
|
||||
for (int i = 1; i < imax + 1; ++i) {
|
||||
P(i, j) += E(i, j);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
Solver copySolver(Solver* solver)
|
||||
{
|
||||
Solver newSolver;
|
||||
newSolver.problem = solver->problem;
|
||||
newSolver.bcLeft = solver->bcLeft;
|
||||
newSolver.bcRight = solver->bcRight;
|
||||
newSolver.bcBottom = solver->bcBottom;
|
||||
newSolver.bcTop = solver->bcTop;
|
||||
newSolver.grid.imax = solver->grid.imax;
|
||||
newSolver.grid.jmax = solver->grid.jmax;
|
||||
newSolver.grid.xlength = solver->grid.xlength;
|
||||
newSolver.grid.ylength = solver->grid.ylength;
|
||||
newSolver.grid.dx = solver->grid.xlength / solver->grid.imax;
|
||||
newSolver.grid.dy = solver->grid.ylength / solver->grid.jmax;
|
||||
newSolver.eps = solver->eps;
|
||||
newSolver.omega = solver->omega;
|
||||
newSolver.itermax = solver->itermax;
|
||||
newSolver.re = solver->re;
|
||||
newSolver.gx = solver->gx;
|
||||
newSolver.gy = solver->gy;
|
||||
newSolver.dt = solver->dt;
|
||||
newSolver.te = solver->te;
|
||||
newSolver.tau = solver->tau;
|
||||
newSolver.gamma = solver->gamma;
|
||||
newSolver.rho = solver->rho;
|
||||
newSolver.levels = solver->levels;
|
||||
newSolver.currentlevel = solver->currentlevel + 1;
|
||||
|
||||
newSolver.r = solver->r;
|
||||
newSolver.e = solver->e;
|
||||
|
||||
return newSolver;
|
||||
}
|
||||
|
||||
void writeResult(Solver* solver)
|
||||
{
|
||||
int imax = solver->grid.imax;
|
||||
int jmax = solver->grid.jmax;
|
||||
double dx = solver->grid.dx;
|
||||
double dy = solver->grid.dy;
|
||||
double* p = solver->p;
|
||||
double* u = solver->u;
|
||||
double* v = solver->v;
|
||||
double x = 0.0, y = 0.0;
|
||||
|
||||
FILE* fp;
|
||||
fp = fopen("pressure.dat", "w");
|
||||
|
||||
if (fp == NULL) {
|
||||
printf("Error!\n");
|
||||
exit(EXIT_FAILURE);
|
||||
}
|
||||
|
||||
for (int j = 1; j < jmax + 1; j++) {
|
||||
y = (double)(j - 0.5) * dy;
|
||||
for (int i = 1; i < imax + 1; i++) {
|
||||
x = (double)(i - 0.5) * dx;
|
||||
fprintf(fp, "%.2f %.2f %f\n", x, y, P(i, j));
|
||||
}
|
||||
fprintf(fp, "\n");
|
||||
}
|
||||
|
||||
fclose(fp);
|
||||
|
||||
fp = fopen("velocity.dat", "w");
|
||||
|
||||
if (fp == NULL) {
|
||||
printf("Error!\n");
|
||||
exit(EXIT_FAILURE);
|
||||
}
|
||||
|
||||
for (int j = 1; j < jmax + 1; j++) {
|
||||
y = dy * (j - 0.5);
|
||||
for (int i = 1; i < imax + 1; i++) {
|
||||
x = dx * (i - 0.5);
|
||||
double vel_u = (U(i, j) + U(i - 1, j)) / 2.0;
|
||||
double vel_v = (V(i, j) + V(i, j - 1)) / 2.0;
|
||||
double len = sqrt((vel_u * vel_u) + (vel_v * vel_v));
|
||||
fprintf(fp, "%.2f %.2f %f %f %f\n", x, y, vel_u, vel_v, len);
|
||||
}
|
||||
}
|
||||
|
||||
fclose(fp);
|
||||
}
|
@ -6,62 +6,21 @@
|
||||
*/
|
||||
#ifndef __SOLVER_H_
|
||||
#define __SOLVER_H_
|
||||
#include "discretization.h"
|
||||
#include "grid.h"
|
||||
#include "parameter.h"
|
||||
|
||||
enum OBJECTBOUNDARY {
|
||||
NONE = 0,
|
||||
TOP,
|
||||
BOTTOM,
|
||||
LEFT,
|
||||
RIGHT,
|
||||
TOPLEFT,
|
||||
BOTTOMLEFT,
|
||||
TOPRIGHT,
|
||||
BOTTOMRIGHT,
|
||||
LOCAL
|
||||
};
|
||||
enum BC { NOSLIP = 1, SLIP, OUTFLOW, PERIODIC };
|
||||
/// @brief
|
||||
enum SHAPE { NOSHAPE = 0, RECT, CIRCLE };
|
||||
|
||||
typedef struct {
|
||||
/* geometry and grid information */
|
||||
Grid grid;
|
||||
/* arrays */
|
||||
double *p, *rhs, **r, **e;
|
||||
double *f, *g;
|
||||
double *u, *v;
|
||||
int* s;
|
||||
Grid* grid;
|
||||
/* parameters */
|
||||
double eps, omega, rho;
|
||||
double re, tau, gamma;
|
||||
double gx, gy;
|
||||
/* time stepping */
|
||||
int itermax, levels, currentlevel;
|
||||
double dt, te;
|
||||
double dtBound;
|
||||
char* problem;
|
||||
int bcLeft, bcRight, bcBottom, bcTop;
|
||||
int itermax;
|
||||
int levels;
|
||||
double **r, **e;
|
||||
} Solver;
|
||||
|
||||
extern void initSolver(Solver*, Parameter*);
|
||||
extern void computeRHS(Solver*);
|
||||
extern double smoothRB(Solver*);
|
||||
extern void restrictMG(Solver*);
|
||||
extern void prolongate(Solver*);
|
||||
extern void correct(Solver*);
|
||||
extern Solver copySolver(Solver*);
|
||||
extern void multiGrid(Solver*);
|
||||
extern void normalizePressure(Solver*);
|
||||
extern void computeTimestep(Solver*);
|
||||
extern void setBoundaryConditions(Solver*);
|
||||
extern void setSpecialBoundaryCondition(Solver*);
|
||||
extern void setObjectBoundaryCondition(Solver*);
|
||||
extern void computeFG(Solver*);
|
||||
extern void adaptUV(Solver*);
|
||||
extern void writeResult(Solver*);
|
||||
extern void print(Solver*, double*);
|
||||
extern void printGrid(Solver*, int*);
|
||||
extern void initSolver(Solver*, Discretization*, Parameter*);
|
||||
extern void solve(Solver*, double*, double*);
|
||||
|
||||
#endif
|
||||
|
@ -7,18 +7,16 @@
|
||||
#include <stdlib.h>
|
||||
#include <time.h>
|
||||
|
||||
double getTimeStamp()
|
||||
double getTimeStamp(void)
|
||||
{
|
||||
struct timespec ts;
|
||||
clock_gettime(CLOCK_MONOTONIC, &ts);
|
||||
return (double)ts.tv_sec + (double)ts.tv_nsec * 1.e-9;
|
||||
}
|
||||
|
||||
double getTimeResolution()
|
||||
double getTimeResolution(void)
|
||||
{
|
||||
struct timespec ts;
|
||||
clock_getres(CLOCK_MONOTONIC, &ts);
|
||||
return (double)ts.tv_sec + (double)ts.tv_nsec * 1.e-9;
|
||||
}
|
||||
|
||||
double getTimeStamp_() { return getTimeStamp(); }
|
||||
|
@ -7,8 +7,7 @@
|
||||
#ifndef __TIMING_H_
|
||||
#define __TIMING_H_
|
||||
|
||||
extern double getTimeStamp();
|
||||
extern double getTimeResolution();
|
||||
extern double getTimeStamp_();
|
||||
extern double getTimeStamp(void);
|
||||
extern double getTimeResolution(void);
|
||||
|
||||
#endif // __TIMING_H_
|
||||
|
@ -7,8 +7,7 @@
|
||||
#ifndef __UTIL_H_
|
||||
#define __UTIL_H_
|
||||
#define HLINE \
|
||||
"------------------------------------------------------------------------" \
|
||||
"----\n"
|
||||
"----------------------------------------------------------------------------\n"
|
||||
|
||||
#ifndef MIN
|
||||
#define MIN(x, y) ((x) < (y) ? (x) : (y))
|
||||
@ -20,4 +19,11 @@
|
||||
#define ABS(a) ((a) >= 0 ? (a) : -(a))
|
||||
#endif
|
||||
|
||||
#define P(i, j) p[(j) * (imax + 2) + (i)]
|
||||
#define F(i, j) f[(j) * (imax + 2) + (i)]
|
||||
#define G(i, j) g[(j) * (imax + 2) + (i)]
|
||||
#define U(i, j) u[(j) * (imax + 2) + (i)]
|
||||
#define V(i, j) v[(j) * (imax + 2) + (i)]
|
||||
#define RHS(i, j) rhs[(j) * (imax + 2) + (i)]
|
||||
|
||||
#endif // __UTIL_H_
|
||||
|
10
EnhancedSolver/2D-seq/vis_files/canal.plot
Normal file
10
EnhancedSolver/2D-seq/vis_files/canal.plot
Normal file
@ -0,0 +1,10 @@
|
||||
unset border; unset tics; unset key;
|
||||
set term gif animate delay 30
|
||||
set output "trace.gif"
|
||||
set xrange [0:30]
|
||||
set yrange [0:4]
|
||||
|
||||
do for [ts=0:120] {
|
||||
plot "particles_".ts.".dat" with points pointtype 7
|
||||
}
|
||||
unset output
|
Loading…
Reference in New Issue
Block a user