---

time_step.c File Reference

#include "utilf.h"
#include <readSIunits.h>
#include "markers.h"
#include "surfaces.h"
#include "extra.h"
#include "stresses.h"
#include "global.h"
#include "printxplot.h"
#include "printfree.h"
#include "momentum.h"
#include "tension.h"
#include "advect.h"
#include "make_bc.h"

Include dependency graph for time_step.c:

Go to the source code of this file.

Functions
real	minau (real2D au, real2D ap, int nx, int ny, int *imin, int *jmin)
real	minav (real2D av, real2D ap, int nx, int ny, int *imin, int *jmin)
void	timestep (int t)

---

Function Documentation

real minau (  real2D    au, real2D    ap, int    nx, int    ny, int *    imin, int *    jmin )

Definition at line 67 of file utilf.c. References INU, real, and real2D. { int i, j; real val = 100.; for (i = 2; i < nx; i++) for (j = 2; j < ny; j++) if (INU(i,j)) if (au[i][j] < val) { val = au[i][j]; *imin = i; *jmin = j; } return val; }

real minav (  real2D    av, real2D    ap, int    nx, int    ny, int *    imin, int *    jmin )

Definition at line 82 of file utilf.c. References INV, real, and real2D. Referenced by timestep(). { int i, j; real val = 100.; for (i = 2; i < nx; i++) for (j = 3; j < ny; j++) if (INV(i,j)) if (av[i][j] < val) { val = av[i][j]; *imin = i; *jmin = j; } return val; }

void timestep (  int    t )

Definition at line 23 of file time_step.c. References advect(), bc_pressure(), bc_vector_bound(), bc_vector_div(), curvature_write(), cut_interface(), divergence(), dp(), findmax(), findmin(), fmodmax(), interface_axisymmetry(), interface_fill(), interface_redis(), interface_rrdz(), interface_rzdr(), interface_splines(), interface_surfaces(), interface_tensionu(), interface_tensionv(), interface_write_markers(), interface_write_splines(), interfaces_axivolume(), interfacial_pressure(), minav(), momentum(), pressure(), printxpl(), rdivmax(), real, reconnect_axis(), sq, sumfield(), sx_fill(), sy_fill(), and viscous_tensor(). { int imax = 0, jmax = 0, iumax = 0, jumax = 0, i1, i, j; real umax, vmax, h, hi2; static real pe, pg; #ifdef SONO static real R, Rp; real a; #endif FILE *msgptr; /*cld*/ real uaxis = 0.0; int iaxis; int test = 0; /*cld*/ h = 1.0 / (num.nx - 2); hi2 = 1.0 / (h * h); /* --------------------------------------------------------------- Control CFL condition ---------------------------------------------------------------- */ if (t % 1 == 0) { umax = fmodmax(u, num.nx, num.ny, &iumax, &jumax); vmax = fmodmax(v, num.nx, num.ny, &i, &j); if (vmax > umax) { umax = vmax; iumax = i; jumax = j; } if (umax > 0.5) { printf("Warning> t=%d velmax=%g x=%d y=%d\n", t, umax, iumax, jumax); printf("CFL number too large. Try smaller tau\n"); printxpl(u, v, ap, p, num.nx, num.ny, num.tau, num.du, num.dp, "cfl.datb"); msgptr = fopen("cfl.xmgr", "wt"); for (i = 0; i < Nint; i++) interface_write_markers(interfaces[i], msgptr); fclose(msgptr); exit(1); } } /* --------------------------------------------------------------- Iteration ---------------------------------------------------------------- */ /*------------------- Advection ------------------- */ for (i = 0; i < Nint; i++) { advect(interfaces[i], u, v, ap, num.nx, num.ny, &uaxis); /* repulse(interfaces[i], 0.5, 1e-5, num.lredis); */ #ifdef RECONNECT_AXIS if (reconnect_axis(&interfaces[i], num.nx, num.ny, 0.1)) { printf("t = %d: reconnect interface %d\n", t, i); msgptr = fopen("reconnect.xmgr", "wt"); interface_write_markers(interfaces[i], msgptr); fclose(msgptr); #if 0 interface_splines(interfaces[i], num.nx, num.ny); printf("Filling fractions ... "); interface_fill(interfaces[i], ap, num.nx, num.ny); for (i1 = 2; i1 <= num.nx; i1++) for (j = 2; j <= num.ny; j++) { au[i1][j] = (ap[i1][j] + ap[i1-1][j])/2.; av[i1][j] = (ap[i1][j] + ap[i1][j-1])/2.; } printf("Ok\n"); fflush(stdout); #endif } #endif interface_splines(interfaces[i], num.nx, num.ny); #ifdef AXISYMMETRIC_SMOOTHING interface_axisymmetry(&interfaces[i], num.lredis, num.nx, num.ny); #else interface_redis(&interfaces[i], num.lredis); #endif interface_splines(interfaces[i], num.nx, num.ny); } pg = phys.po*exp(phys.gamma*log(phys.vo/ fabs(interfaces_axivolume(interfaces, Nint)))); pe = phys.pext - phys.pe*sin(phys.omega*phys.t); a = exp(1/3.*log(fabs(interfaces_axivolume(interfaces, Nint))*3./(4.*PI))); if (t == 0) Rp = 0.0; else Rp = a - R; R = a; /* compute interfacial pressure */ interfacial_pressure(interfaces[0], u, v, ap, phys.sigma, phys.visliq, pg, num.tau, num.nx, num.ny); cut_interface(interfaces[0], &cutp, 0.0, 0.0); cut_interface(interfaces[0], &cutu, -0.5, 0.0); cut_interface(interfaces[0], &cutv, 0.0, -0.5); interface_surfaces(interfaces[0], cutp, sxp, syp, ap, cc, S12, work, num.nx, num.ny); interface_surfaces(interfaces[0], cutu, sxu, syu, au, a1, S12, work, num.nx, num.ny); interface_surfaces(interfaces[0], cutv, sxv, syv, av, a2, S12, work, num.nx, num.ny); sx_fill(sxp, ap, 0.0, 0.0, interfaces[0], num.nx, num.ny); sy_fill(syp, ap, 0.0, 0.0, interfaces[0], num.nx, num.ny); sx_fill(sxu, au, -0.5, 0.0, interfaces[0], num.nx, num.ny); sy_fill(syu, au, -0.5, 0.0, interfaces[0], num.nx, num.ny); for (j = 2; j <= num.ny; j++) syu[num.nx+1][j] = (real)j - 1.5; sx_fill(sxv, av, 0.0, -0.5, interfaces[0], num.nx, num.ny); for (i = 2; i <= num.nx; i++) sxv[i][num.ny+1] = (real)num.ny - 1.5; sy_fill(syv, av, 0.0, -0.5, interfaces[0], num.nx, num.ny); bc_vector_div(u, v, ap, interfaces[0], num.nx, num.ny); /*--------------------- Momentum equation -------------------------*/ viscous_tensor(u, v, S11, S22, S12, phys.visliq, num.tau, num.nx, num.ny); momentum(u, v, a1, a2, ap, av, sxu, syu, sxv, syv, interfaces[0], cutu, cutv, rz1dr, rz2dr, rr1dz, rr2dz, work, S11, S22, S12, divs, phys.visliq, phys.gx, phys.gy, num.tau, num.nx, num.ny); interface_tensionu(interfaces[0], cutu, u, v, a1, ap, syu, p, S11, S22, num.nx, num.ny); interface_tensionv(interfaces[0], cutv, u, v, a2, ap, sxv, av, p, S11, S22, num.nx, num.ny); bc_vector_bound(u, v, num.nx, num.ny); bc_pressure(p, ap, pe*sq(num.tau/h), R, num.nx, num.ny); interface_rzdr(interfaces[0], cutp, -0.5, rz1dr, num.nx, num.ny); interface_rzdr(interfaces[0], cutp, 0.5, rz2dr, num.nx, num.ny); interface_rrdz(interfaces[0], cutp, -0.5, rr1dz, num.nx, num.ny); interface_rrdz(interfaces[0], cutp, 0.5, rr2dz, num.nx, num.ny); pressure(u, v, p, a1, a2, ap, av, sxp, syp, syu, sxv, rz1dr, rz2dr, rr1dz, rr2dz, divs, work1, S11, S22, S12, a3, work, interfaces[0], num.lredis, &num.cycles, num.npre, num.npost, num.maxdiv, num.nx, num.ny, num.ng); bc_vector_div(u, v, ap, interfaces[0], num.nx, num.ny); /* --------------------------------------------------------------- Outputs ---------------------------------------------------------------- */ if (t % out.txplot == 0) { char s[256]; /* sprintf(s, "free.%d.fig", t / out.txplot); msgptr = fopen(s, "wt"); printfree(msgptr, u, v, ap, interfaces[0], num.nx, num.ny); fclose(msgptr); */ sprintf(s, "data.%d.datb", t/out.txplot); printxpl(u, v, ap, p, num.nx, num.ny, num.tau, num.du, num.dp, s); } if (t % out.tprint == 0) { if (t == 0) msgptr = fopen("msg", "wt"); else msgptr = fopen("msg", "at"); printf("t=%d umax=%g at (%d,%d) ncycle=%d", t, umax, iumax, jumax, num.cycles); if (phys.sigma > 0.0) printf(" cwn=%g\n", num.tau/sqrt(h*h*h/phys.sigma)); else printf("\n"); divergence(u, v, ap, sxp, syp, rz1dr, rz2dr, rr1dz, rr2dz, work, num.nx, num.ny); printf("t=%d div max before=%g div max after=%g\n", t, rdivmax(divs, num.nx, num.ny, &imax, &jmax), rdivmax(work, num.nx, num.ny, &imax, &jmax)); fprintf(msgptr, "%d %g %g %g %g %g %g %g %g %g %g", t, phys.t, umax * h / num.tau * num.du, interfaces_axivolume(interfaces, Nint), findmax(ap, num.nx, num.ny, &imax, &jmax), rdivmax(work, num.nx, num.ny, &imax, &jmax), num.tau, interfaces_x_amp(interfaces, Nint)/(num.nx - 2), interfaces_y_amp(interfaces, Nint)/(num.ny - 2), /* interfaces_ymoment(interfaces, Nint), */ pg*num.dp, pe*num.dp); if (interfaces[0].n - 1 % 2 == 0) { i = (interfaces[0].n - 1)/2; fprintf(msgptr, " %g %g\n", sqrt(sq(interfaces[0].x[i] - 0.5*(num.nx - 2) - 2.0) + sq(interfaces[0].y[i] - 2.0)) - exp(1./3.*log(3./(4.*PI)*fabs(interfaces_axivolume(interfaces, Nint)))), exp(1./3.*log(3./(4.*PI)*fabs(interfaces_axivolume(interfaces, Nint))))); } else { i = (interfaces[0].n - 1)/2; fprintf(msgptr, " %g %g\n", 0.5*(sqrt(sq(interfaces[0].x[i] - 0.5*(num.nx - 2) - 2.0) + sq(interfaces[0].y[i] - 2.0)) + sqrt(sq(interfaces[0].x[i+1] - 0.5*(num.nx - 2) - 2.0) + sq(interfaces[0].y[i+1] - 2.0))) - exp(1./3.*log(3./(4.*PI)*fabs(interfaces_axivolume(interfaces, Nint)))), exp(1./3.*log(3./(4.*PI)*fabs(interfaces_axivolume(interfaces, Nint))))); } printf("mass=%g\n", sumfield(cc, num.nx, num.ny)); printf("min c: %g max c: %g\n", findmin(ap, num.nx, num.ny, &imax, &jmax), findmax(ap, num.nx, num.ny, &imax, &jmax)); a = minau(au, ap, num.nx, num.ny, &imax, &jmax); printf("min au: %f at %d,%d\n", a, imax, jmax); a = minav(av, ap, num.nx, num.ny, &imax, &jmax); printf("min av: %f at %d,%d\n", a, imax, jmax); printf("dp (SI): %g\n", dp(p, cc, num.nx, num.ny) * h * h / sq(num.tau) * num.dp); fflush(stdout); fclose(msgptr); } if (t % out.tmark == 0) { if (t == 0) msgptr = fopen("markers.xmgr", "wt"); else msgptr = fopen("markers.xmgr", "at"); printf("volume: %g Nint: %d\n", interfaces_axivolume(interfaces, Nint), Nint); for (i = 0; i < Nint; i++) interface_write_markers(interfaces[i], msgptr); // interface_write_markers_speed(phys.t, interfaces[i], u, v, // num.tau, num.du, num.nx, num.ny, msgptr); fclose(msgptr); #if 0 if (t == 0) msgptr = fopen("splines.xmgr", "wt"); else msgptr = fopen("splines.xmgr", "at"); for (i = 0; i < Nint; i++) interface_write_splines(interfaces[i], msgptr, 1000); fclose(msgptr); #endif if (t == 0) msgptr = fopen("curvature.xmgr", "wt"); else msgptr = fopen("curvature.xmgr", "at"); for (i = 0; i < Nint; i++) curvature_write(interfaces[i], msgptr); fclose(msgptr); if (t == 0) msgptr = fopen("ipressure.xmgr", "wt"); else msgptr = fopen("ipressure.xmgr", "at"); for (i = 0; i < interfaces[0].n; i++) fprintf(msgptr, "%g %g\n", interfaces[0].t[i]/interfaces[0].t[interfaces[0].n-1], interfaces[0].p[i]); fprintf(msgptr, "\n"); fclose(msgptr); if (t == 0) msgptr = fopen("radius.xmgr", "wt"); else msgptr = fopen("radius.xmgr", "at"); for (i = 0; i < Nint; i++) { fprintf(msgptr, "\n"); for (j = 0; j < interfaces[i].n - 1; j++) fprintf(msgptr, "%g %g\n", interfaces[i].t[j]/ interfaces[i].t[interfaces[i].n - 1], sqrt(sq(interfaces[i].x[j] - 0.5*(num.nx - 2) - 2.0) + sq(interfaces[i].y[j] - 2.0))); } fclose(msgptr); /*cld*/ if (t == 0) msgptr = fopen("axis.xmgr", "wt"); else msgptr = fopen("axis.xmgr", "at"); iaxis = interfaces[0].x[0]; // printf("%g %d \n", phys.t*num.tausi, iaxis); fprintf(msgptr, "%g %g %g %g %g\n", phys.t, (interfaces[0].x[0]-2.0)*h, u[iaxis-1][1] * h / num.tau * num.du, u[iaxis-1][2] * h / num.tau * num.du, u[iaxis-1][3] * h / num.tau * num.du); fclose(msgptr); if (t == 0) msgptr = fopen("velocity.xmgr", "wt"); else msgptr = fopen("velocity.xmgr", "at"); // if ((interfaces[0].x[0]-2.0)*h > 0.616406 && test == 0){ /* Double Size Bubble */ // if ((interfaces[0].x[0]-2.0)*h > 0.918281 && test == 0){ /* Double Size Bubble */ iaxis = interfaces[0].x[0]; test = 1; fprintf(msgptr, "%g %g %g %g %g\n", phys.t, (interfaces[0].x[0]-2.0)*h, u[iaxis-1][1] * h / num.tau * num.du, u[iaxis-1][2] * h / num.tau * num.du, u[iaxis-1][3] * h / num.tau * num.du); // exit(1); // } fclose(msgptr); /*cld*/ } phys.t += num.tau; /*cld for (j = 1; j < interfaces[0].n ; j++){ if (interfaces[0].y[j] <= 3.0 && interfaces[0].y[j] < interfaces[0].y[j-1] && interfaces[0].x[j] < interfaces[0].x[j-1]) exit(1); } cld*/ #ifndef NO_ADAPTATIVE adaptative(u, v, p, &num.tau, num.cycles, num.nx, num.ny); #endif }

---