---

momentum.c File Reference

#include "utilf.h"
#include "markers.h"
#include "surfaces.h"
#include "extra.h"
#include "momentum.h"
#include "mg.h"
#include "make_bc.h"
#include "fluxes.h"
#include "inout.h"

Include dependency graph for momentum.c:

Go to the source code of this file.

Functions
void	momentum (real2D u, real2D v, real2D cu, real2D cv, real2D ap, real2D av, real2D sxu, real2D syu, real2D sxv, real2D syv, interface in, interface_cut cutu, interface_cut cutv, real2D rz1dr, real2D rz2dr, real2D rr1dz, real2D rr2dz, real2D Sx12, real2D S11, real2D S22, real2D S12, real2D Sy12, real visliq, real gx, real gy, real tau, int nx, int ny)
void	pressure (real2D u, real2D v, real2D p, real2D cu, real2D cv, real2D ap, real2D av, real2D sxp, real2D syp, real2D syu, real2D sxv, real2D rz1dr, real2D rz2dr, real2D rr1dz, real2D rr2dz, real2D div, real2D tmp, real2D c0, real2D c1, real2D c2, real2D c3, real2D c4, interface in, real lredis, int *ncycle, int npre, int npost, real maxdiv, int nx, int ny, int ng)
void	coli (real2D sxp, real2D syp, real2D syu, real2D sxv, real2D cu, real2D cv, real2D ap, real2D av, real2D rz1dr, real2D rz2dr, real2D rr1dz, real2D rr2dz, real2D c0, real2D c1, real2D c2, real2D c3, real2D c4, int nx, int ny, real scale)
void	adaptative (real2D u, real2D v, real2D p, real *tau, int ncycle, int nx, int ny)
void	rescale (real2D u, real2D v, real2D p, real tscale, int nx, int ny)
void	divergence (real2D u, real2D v, real2D ap, real2D sxp, real2D syp, real2D rz1dr, real2D rz2dr, real2D rr1dz, real2D rr2dz, real2D div, int nx, int ny)

---

Function Documentation

void adaptative (  real2D    u, real2D    v, real2D    p, real *    tau, int    ncycle, int    nx, int    ny )

Definition at line 321 of file momentum.c. { real umax, vmax; int i, j; umax = fmodmax(u, nx, ny, &i, &j); vmax = fmodmax(v, nx, ny, &i, &j); if (vmax > umax) umax = vmax; if (ncycle > MAXCYCLE || umax > MAXVEL) { rescale(u, v, p, 1./DIVIDE, nx, ny); *tau /= DIVIDE; } else if (ncycle < 2 && umax*MULT < MAXVEL) { rescale(u, v, p, MULT, nx, ny); *tau *= MULT; } }

void coli (  real2D    sxp, real2D    syp, real2D    syu, real2D    sxv, real2D    cu, real2D    cv, real2D    ap, real2D    av, real2D    rz1dr, real2D    rz2dr, real2D    rr1dz, real2D    rr2dz, real2D    c0, real2D    c1, real2D    c2, real2D    c3, real2D    c4, int    nx, int    ny, real    scale )

Definition at line 252 of file momentum.c. References INP, real, and real2D. Referenced by mgsolve(), and pressure(). { int i, j; real syp1, syp3, sxp2, sxp4; for (i = 2; i < nx; i++) for (j = 2; j < ny; j++) if (INP(i,j)) { syp1 = scale*(syp[i+1][j] + rz1dr[i][j])/cu[i+1][j]; syp3 = scale*(syp[i][j] + rz2dr[i][j])/cu[i][j]; c1[i][j] = syu[i+2][j]*syp1; c3[i][j] = syu[i][j]*syp3; sxp2 = j == ny-1 ? 0.0 : scale*(sxp[i][j+1] - rr1dz[i][j])/cv[i][j+1]; sxp4 = j == 2 ? 0.0 : scale*(sxp[i][j] - rr2dz[i][j])/cv[i][j]; c2[i][j] = (sxv[i][j+2] - 0.5*av[i][j+1])*sxp2; c4[i][j] = (sxv[i][j] + 0.5*av[i][j])*sxp4; c0[i][j] = syu[i+1][j]*(syp1 + syp3) + sxv[i][j+1]*(sxp2 + sxp4) + 0.5*(av[i][j+1]*sxp2 - av[i][j]*sxp4); #ifdef CHECK_DOMINANCE if (fabs(c0[i][j]) < fabs(c1[i][j])*INP(i+1,j) + fabs(c2[i][j])*INP(i,j+1) + fabs(c3[i][j])*INP(i-1,j) + fabs(c4[i][j])*INP(i,j-1)) { printf("WARNING: coli(%d,%d): c0: %g c1: %g c2: %g c3: %g c4: %g\n", i, j, c0[i][j], c1[i][j]*INP(i+1,j), c2[i][j]*INP(i,j+1), c3[i][j]*INP(i-1,j), c4[i][j]*INP(i,j-1)); c0[i][j] = fabs(c1[i][j]) + fabs(c2[i][j]) + fabs(c3[i][j]) + fabs(c4[i][j]); printf(" setting c0 to %g\n", c0[i][j]); } #endif #if 0 if (i == 150 && j == 3 && scale == 1.0) { printf("syp[%d][%d]: %g rz1dr[%d][%d]: %g cu[%d][%d]: %g\n", i+1, j, syp[i+1][j], i, j, rz1dr[i][j], i+1, j, cu[i+1][j]); printf("syp[%d][%d]: %g rz2dr[%d][%d]: %g cu[%d][%d]: %g\n", i, j, syp[i][j], i, j, rz2dr[i][j], i, j, cu[i][j]); printf("syu[%d][%d]: %g c1[%d][%d]: %g\n", i+2, j, syu[i+2][j], i, j, c1[i][j]); printf("syu[%d][%d]: %g c3[%d][%d]: %g\n", i, j, syu[i][j], i, j, c3[i][j]); printf("\n"); printf("sxp[%d][%d]: %g rr1dz[%d][%d]: %g cv[%d][%d]: %g\n", i, j+1, sxp[i][j+1], i, j, rr1dz[i][j], i, j+1, cv[i][j+1]); printf("sxp[%d][%d]: %g rr2dz[%d][%d]: %g cv[%d][%d]: %g\n", i, j, sxp[i][j], i, j, rr2dz[i][j], i, j, cv[i][j]); printf("sxv[%d][%d]: %g av[%d][%d]: %g c2[%d][%d]: %g\n", i, j+2, sxv[i][j+2], i, j+1, av[i][j+1], i, j, c2[i][j]); printf("sxv[%d][%d]: %g av[%d][%d]: %g c4[%d][%d]: %g\n", i, j, sxv[i][j], i, j, av[i][j], i, j, c4[i][j]); printf("\n"); printf("syu[%d][%d]: %g sxv[%d][%d]: %g c0[%d][%d]: %g\n", i+1, j, syu[i+1][j], i, j+1, sxv[i][j+1], i, j, c0[i][j]); } #endif } }

void divergence (  real2D    u, real2D    v, real2D    ap, real2D    sxp, real2D    syp, real2D    rz1dr, real2D    rz2dr, real2D    rr1dz, real2D    rr2dz, real2D    div, int    nx, int    ny )

Definition at line 358 of file momentum.c. References INP, and real2D. Referenced by pressure(), and timestep(). { int i, j; for (i = 2; i < nx; i++) for (j = 2; j < ny; j++) if (INP(i,j)) { div[i][j] = (syp[i+1][j] + rz1dr[i][j])*u[i+1][j] - (syp[i][j] + rz2dr[i][j])*u[i][j] + (sxp[i][j+1] - rr1dz[i][j])*v[i][j+1] - (sxp[i][j] - rr2dz[i][j])*v[i][j]; } else div[i][j] = 0.0; }

void momentum (  real2D    u, real2D    v, real2D    cu, real2D    cv, real2D    ap, real2D    av, real2D    sxu, real2D    syu, real2D    sxv, real2D    syv, interface    in, interface_cut    cutu, interface_cut    cutv, real2D    rz1dr, real2D    rz2dr, real2D    rr1dz, real2D    rr2dz, real2D    Sx12, real2D    S11, real2D    S22, real2D    S12, real2D    Sy12, real    visliq, real    gx, real    gy, real    tau, int    nx, int    ny )

Definition at line 65 of file momentum.c. { int i, j; real h, f; h = 1./(nx - 2); gx *= sq(tau)/h; gy *= sq(tau)/h; visliq *= tau/(h*h); /* ----------------------------------------------------------------- Advection term ------------------------------------------------------------------ */ #ifdef UPWIND for (i = 1; i <= nx; i++) for (j = 2; j < ny; j++) { f = 0.5*(u[i][j] + u[i+1][j]); S11[i][j] -= f > 0.0 ? u[i][j]*f : u[i+1][j]*f; } for (i = 2; i < nx; i++) for (j = 2; j <= ny; j++) { f = 0.5*(v[i][j] + v[i][j+1]); S22[i][j] -= f > 0.0 ? v[i][j]*f : v[i][j+1]*f; } for (i = 2; i <= nx; i++) for (j = 2; j <= ny; j++) { f = 0.5*(v[i][j] + v[i-1][j]); Sx12[i][j] = S12[i][j] - (f > 0.0 ? u[i][j-1]*f : u[i][j]*f); f = 0.5*(u[i][j] + u[i][j-1]); Sy12[i][j] = S12[i][j] - (f > 0.0 ? v[i-1][j]*f : v[i][j]*f); } #else for (i = 1; i <= nx; i++) for (j = 2; j < ny; j++) S11[i][j] -= 0.25*sq(u[i][j] + u[i+1][j]); for (i = 2; i < nx; i++) for (j = 2; j <= ny; j++) S22[i][j] -= 0.25*sq(v[i][j] + v[i][j+1]); for (i = 2; i <= nx; i++) for (j = 2; j <= ny; j++) Sx12[i][j] = Sy12[i][j] = S12[i][j] - 0.25*(u[i][j] + u[i][j-1])*(v[i][j] + v[i-1][j]); #endif /* ------------------------------------------------------------------ div(S) Gravity term ------------------------------------------------------------------ */ interface_rzdr(in, cutu, -0.5, rz1dr, nx, ny); interface_rzdr(in, cutu, 0.5, rz2dr, nx, ny); interface_rrdz(in, cutu, -0.5, rr1dz, nx, ny); interface_rrdz(in, cutu, 0.5, rr2dz, nx, ny); for (i = 2; i <= nx; i++) for (j = 2; j < ny; j++) if (INU(i,j)) u[i][j] += ((syu[i+1][j] + rz1dr[i][j])*S11[i][j] - (syu[i][j] + rz2dr[i][j])*S11[i-1][j] + (sxu[i][j+1] - rr1dz[i][j])*Sx12[i][j+1] - (sxu[i][j] - rr2dz[i][j])*Sx12[i][j])/cu[i][j] + gx; interface_rzdr(in, cutv, -0.5, rz1dr, nx, ny); interface_rzdr(in, cutv, 0.5, rz2dr, nx, ny); interface_rrdz(in, cutv, -0.5, rr1dz, nx, ny); interface_rrdz(in, cutv, 0.5, rr2dz, nx, ny); for (i = 2; i < nx; i++) for (j = 3; j <= ny; j++) if (INV(i,j)) v[i][j] += ((syv[i+1][j] + rz1dr[i][j])*Sy12[i+1][j] - (syv[i][j] + rz2dr[i][j])*Sy12[i][j] + (sxv[i][j+1] - rr1dz[i][j])*S22[i][j] - (sxv[i][j] - rr2dz[i][j])*S22[i][j-1] - av[i][j]*2.*visliq*v[i][j]/((real)j - 2.))/cv[i][j] + gy; }

void pressure (  real2D    u, real2D    v, real2D    p, real2D    cu, real2D    cv, real2D    ap, real2D    av, real2D    sxp, real2D    syp, real2D    syu, real2D    sxv, real2D    rz1dr, real2D    rz2dr, real2D    rr1dz, real2D    rr2dz, real2D    div, real2D    tmp, real2D    c0, real2D    c1, real2D    c2, real2D    c3, real2D    c4, interface    in, real    lredis, int *    ncycle, int    npre, int    npost, real    maxdiv, int    nx, int    ny, int    ng )

Definition at line 153 of file momentum.c. { int i, j; /* ------------------------------------------------------------------ Resolution of the pressure equation div=divergence p=pressure ------------------------------------------------------------------ */ divergence(u, v, ap, sxp, syp, rz1dr, rz2dr, rr1dz, rr2dz, div, nx, ny); #if 0 datbarray(div, nx, ny, "divbefore"); datbarray(u, nx, ny, "ubefore"); datbarray(v, nx, ny, "vbefore"); datbarray(p, nx, ny, "pbefore"); datbarray(rz1dr, nx, ny, "rz1dr"); datbarray(rz2dr, nx, ny, "rz2dr"); datbarray(rr1dz, nx, ny, "rr1dz"); datbarray(rr2dz, nx, ny, "rr2dz"); datbarray(sxp, nx, ny, "sxp"); datbarray(syp, nx, ny, "syp"); datbarray(sxv, nx, ny, "sxv"); datbarray(syu, nx, ny, "syu"); datbarray(ap, nx, ny, "ap"); datbarray(cv, nx, ny, "cv"); datbarray(av, nx, ny, "av"); datbarray(cu, nx, ny, "cu"); { FILE *fptr = fopen("sxp", "wt"); sx_write(sxp, 0.0, 0.0, nx, ny, fptr); fclose(fptr); fptr = fopen("sxv", "wt"); sx_write(sxv, 0.0, -0.5, nx, ny, fptr); fclose(fptr); fptr = fopen("toto", "wt"); interface_write_markers(in, fptr); fclose(fptr); } fill0(c0, nx, ny); fill0(c1, nx, ny); fill0(c2, nx, ny); fill0(c3, nx, ny); fill0(c4, nx, ny); #endif coli(sxp, syp, syu, sxv, cu, cv, ap, av, rz1dr, rz2dr, rr1dz, rr2dz, c0, c1, c2, c3, c4, nx, ny, 1.0); #if 0 datbarray(c0, nx, ny, "c0"); datbarray(c1, nx, ny, "c1"); datbarray(c2, nx, ny, "c2"); datbarray(c3, nx, ny, "c3"); datbarray(c4, nx, ny, "c4"); #endif #if 1 copy(tmp, p, nx, ny); *ncycle = mgsolve(div, p, ap, c0, c1, c2, c3, c4, in, lredis, maxdiv, nx, ny, ng); #else for (i = 0; i < 20; i++) { relax(p, div, ap, c0, c1, c2, c3, c4, nx, ny); bc_scalar(p, nx, ny, NULGRAD2); } #endif /* ------------------------------------------------------------------ Pressure correction ------------------------------------------------------------------ */ for (i = 2; i <= nx; i++) for (j = 2; j < ny; j++) if (INU(i,j)) u[i][j] += (syu[i][j]*p[i-1][j] - syu[i+1][j]*p[i][j])/cu[i][j]; for (i = 2; i < nx; i++) for (j = 3; j < ny; j++) /* !!!! <= ny when No top wall !!!! */ if (INV(i,j)) v[i][j] += ((sxv[i][j] + av[i][j]/2.)*p[i][j-1] - (sxv[i][j+1] - av[i][j]/2.)*p[i][j])/cv[i][j]; #if 0 datbarray(u, nx, ny, "uafter"); datbarray(v, nx, ny, "vafter"); datbarray(p, nx, ny, "pafter"); #endif }

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

Definition at line 342 of file momentum.c. { int i, j; real tscale2 = sq(tscale); for (i = 1; i <= nx; i++) for (j = 1; j <= ny; j++) { if (u[i][j] != UNDEFINED) u[i][j] *= tscale; if (v[i][j] != UNDEFINED) v[i][j] *= tscale; p[i][j] *= tscale2; } }

---