---

momentum.denis.c File Reference

#include "utilf.h"
#include "markers.h"
#include "extra.h"
#include "momentum.h"
#include "vmg.h"
#include "make_bc.h"

Include dependency graph for momentum.denis.c:

Go to the source code of this file.

Functions
void	momentum (real2D u, real2D v, real2D a, real2D c, real2D w1, real2D S11, real2D S22, real2D S12, real2D w2, real visliq, real gx, real gy, real tau, int nx, int ny)
void	pressure (real2D u, real2D v, real2D p, real2D cc, real2D a, real2D c, real2D div, real2D nu1, real2D nu2, real2D nu3, real2D nu4, real2D c0, real2D c1, real2D c2, real2D c3, real2D c4, real epsilon, int nx, int ny, int ng)
void	adaptative (real2D u, real2D v, real2D divafter, real2D cc, int *npre, int *npost, real mineps, real maxeps, int nx, int ny)
void	adaptative_tau (real2D u, real2D v, real2D divafter, real2D cc, real2D p, real *tau, real mineps, real maxeps, int nx, int ny)
void	rescale (real2D u, real2D v, real2D p, real tscale, int nx, int ny)

---

Function Documentation

void adaptative (  real2D    u, real2D    v, real2D    divafter, real2D    cc, int *    npre, int *    npost, real    mineps, real    maxeps, int    nx, int    ny )

Definition at line 281 of file momentum.denis.c. References MAX, rdivmax(), real, and real2D. { int i, j; real epsilon, h = 1.0/(nx - 2), hi2; hi2 = 1.0/(h*h); /* ------------------------------------------------------------------ Adaptation of the multigrid for the given precision ------------------------------------------------------------------ */ for (i = 2; i <= nx - 1; i++) for (j = 2; j <= ny - 1; j++) if (IN(cc[i][j])) divafter[i][j] = (((real)j - 1.)*v[i][j+1] - ((real)j - 2.)*v[i][j] + ((real)j - 1.5)*(u[i+1][j] - u[i][j]))*hi2; else divafter[i][j] = 0.0; /* datbarray(divafter, nx, ny, "divafter"); */ epsilon = rdivmax(divafter, nx, ny); if (epsilon > mineps) { (*npre)++; (*npost)++; } else if (epsilon < maxeps) { *npre = MAX(*npre - 1, 1); *npost = MAX(*npost - 1, 1); } }

void adaptative_tau (  real2D    u, real2D    v, real2D    divafter, real2D    cc, real2D    p, real *    tau, real    mineps, real    maxeps, int    nx, int    ny )

Definition at line 313 of file momentum.denis.c. References datbarray(), findninf(), rdivmax(), real, real2D, and rescale(). { int i, j; real epsilon, velmax, h = 1.0/(nx - 2), hi2; hi2 = 1./(h*h); /* ------------------------------------------------------------------ Adaptation of the time step for the given precision ------------------------------------------------------------------ */ for (i = 2; i <= nx - 1; i++) for (j = 2; j <= ny - 1; j++) if (IN(cc[i][j])) divafter[i][j] = (((real)j - 1.)*v[i][j+1] - ((real)j - 2.)*v[i][j] + ((real)j - 1.5)*(u[i+1][j] - u[i][j]))*hi2; else divafter[i][j] = 0.0; datbarray(divafter, nx, ny, "divafter"); epsilon = rdivmax(divafter, nx, ny); velmax = findninf(u, v, nx, ny, &i, &j); if (epsilon > mineps || velmax > 0.1) { rescale(u, v, p, 0.5, nx, ny); *tau *= 0.5; } else if (epsilon < maxeps && velmax < 0.05) { rescale(u, v, p, 2.0, nx, ny); *tau *= 2.; } }

void momentum (  real2D    u, real2D    v, real2D    a, real2D    c, real2D    w1, real2D    S11, real2D    S22, real2D    S12, real2D    w2, real    visliq, real    gx, real    gy, real    tau, int    nx, int    ny )

Definition at line 70 of file momentum.denis.c. References real, real2D, and sq. { int i, j; real h; real gx2, gy2, tau2h; h = 1. / (nx - 2); tau2h = sq(tau) / h; gx2 = gx * tau2h; gy2 = gy * tau2h; visliq *= tau/(h*h); /* ----------------------------------------------------------------- Advection term ------------------------------------------------------------------ */ for (i = 2; i <= nx - 1; i++) { for (j = 2; j <= ny - 1; j++) if (IN(a[i-1][j])) w1[i][j] = 0.25*( (u[i][j] + u[i][j-1])*(v[i][j] + v[i-1][j]) - (u[i][j+1] + u[i][j])*(v[i][j+1] + v[i-1][j+1]) + sq(u[i-1][j] + u[i][j]) - sq(u[i][j] + u[i+1][j]) /* axisymetric term */ - (v[i][j] + v[i][j+1] + v[i-1][j+1] + v[i-1][j]) *u[i][j]/((real)j - 1.5)); for (j = 3; j <= ny - 1; j++) if (IN(c[i][j-1])) w2[i][j] = 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+1][j] + v[i][j]) + sq(v[i][j-1] + v[i][j]) - sq(v[i][j] + v[i][j+1]) /* axisymetric term */ - sq(2.0*v[i][j])/((real)j - 2.0)); } /* ------------------------------------------------------------------ div(S) Gravity term ------------------------------------------------------------------ */ for (i = 2; i <= nx - 1; i++) { for (j = 2; j <= ny - 1; j++) if (IN(a[i-1][j])) u[i][j] += w1[i][j] + S12[i][j+1] - S12[i][j] + S11[i][j] - S11[i-1][j] /* axisymetric term */ /* + 0.5*(S12[i][j] + S12[i][j+1])/((real)j - 1.5) */ + gx2; for (j = 3; j <= ny - 1; j++) if (IN(c[i][j-1])) v[i][j] += w2[i][j] + S12[i+1][j] - S12[i][j] + S22[i][j] - S22[i][j-1] /* axisymetric term */ + (0.5*(S22[i][j] + S22[i][j-1]) - 2.*visliq*v[i][j]/ ((real)j - 2.0)) /((real)j - 2.0) + gy2; } }

void pressure (  real2D    u, real2D    v, real2D    p, real2D    cc, real2D    a, real2D    c, real2D    div, real2D    nu1, real2D    nu2, real2D    nu3, real2D    nu4, real2D    c0, real2D    c1, real2D    c2, real2D    c3, real2D    c4, real    epsilon, int    nx, int    ny, int    ng )

Definition at line 143 of file momentum.denis.c. References bc_scalar(), coli(), datbarray(), mglin(), NULGRAD2, real, real2D, and relax(). { int i, j; real h = 1.0/(nx - 2), hi2; real p0, p3, p4, n3, n4; hi2 = 1./(h*h); /* ------------------------------------------------------------------ Resolution of the pressure equation div=divergence p=pressure ------------------------------------------------------------------ */ for (i = 2; i <= nx - 1; i++) for (j = 2; j <= ny - 1; j++) if (IN(cc[i][j])) div[i][j] = (((real)j - 1.)*v[i][j+1] - ((real)j - 2.)*v[i][j] + ((real)j - 1.5)*(u[i+1][j] - u[i][j])) #ifndef MG *hi2 #endif ; else div[i][j] = 0.0; datbarray(div, nx, ny, "divbefore"); datbarray(p, nx, ny, "pbefore"); datbarray(cc, nx, ny, "cc"); datbarray(u, nx, ny, "ubefore"); datbarray(v, nx, ny, "vbefore"); coli(cc, nu1, nu2, nu3, nu4, c0, c1, c2, c3, c4, nx, ny); #ifdef MG /* mginit(div, p, cc, ra, rc, rei, pint, work, in, nx, ny, ng); for (i = 0; i < ncycle; i++) mgvcycle(ng, nx, ny, npre, npost); */ for (i = 0; i < ncycle; i++) { relax(p, div, cc, nu1, nu2, nu3, nu4, c0, c1, c2, c3, c4, nx, ny, pa); bc_scalar(p, nx, ny, NULGRAD2); } #else mglin(div, p, cc, ra, rc, rei, epsilon, nx, ny, ng); #endif datbarray(p, nx, ny, "pafter"); /* ------------------------------------------------------------------ Pressure correction ------------------------------------------------------------------ */ for (i = 2; i <= nx - 1; i++) { for (j = 2; j <= ny - 1; j++) if (IN(a[i-1][j])) { if (IN(cc[i][j])) { p0 = p[i][j]; if (IN(cc[i-1][j])) { p3 = p[i-1][j]; n3 = 1.0; } else { p3 = pa; n3 = -nu3[i][j]; } } else if (IN(cc[i-1][j])) { p3 = p[i-1][j]; p0 = pa; n3 = nu1[i-1][j]; } else { fprintf(stderr, "pressure(): error for u(%d,%d)\n", i, j); exit(1); } u[i][j] -= (p0 - p3)/n3; } for (j = 3; j <= ny - 1; j++) if (IN(c[i][j-1])) { if (IN(cc[i][j])) { p0 = p[i][j]; if (IN(cc[i][j-1])) { p4 = p[i][j-1]; n4 = 1.0; } else { p4 = pa; n4 = -nu4[i][j]; } } else if (IN(cc[i][j-1])) { p4 = p[i][j-1]; p0 = pa; n4 = nu2[i][j-1]; } else { fprintf(stderr, "pressure(): error for v(%d,%d)\n", i, j); exit(1); } v[i][j] -= (p0 - p4)/n4; } } }

void rescale (  real2D    u, real2D    v, real2D    p, real    tscale, int    nx, int    ny )

Definition at line 348 of file momentum.denis.c. References real, real2D, sq, and UNDEFINED. Referenced by adaptative(), and adaptative_tau(). { int i, j; real tscale2 = sq(tscale); for (i = 1; i <= nx; i++) for (j = 1; j <= ny; j++) { u[i][j] *= tscale; v[i][j] *= tscale; p[i][j] *= tscale2; } }

---