---

time_step.simple.c

Go to the documentation of this file.

/*------------------------------------------------------------------*/
/* SCCS Information: %W%    %G% */
/*------------------------------------------------------------------*/
#include "utilf.h"
#include <readSIunits.h>
#include "markers.h"
#include "extra.h"
#include "stresses.h"
#include "global.h"
#include "printxplot.h"
#include "momentum.h"
#include "make_bc.h"

FILE *interu, *interv;

void timestep(int t)
{
  int imax = 0, jmax = 0, i, j;
  real velmax, h, hi2;
  static real pe, pext, pg;
#ifdef SONO
  static real R, Rp, pR;
  real a;
#endif
  real diameter;
  FILE *msgptr;

  h = 1.0 / (num.nx - 2);
  hi2 = 1.0 / (h * h);
  
  /* ---------------------------------------------------------------
     Control CFL condition
     ---------------------------------------------------------------- */
#if 0
  if(t % 1 == 0) {
    velmax = findninf(u, v, num.nx, num.ny, &imax, &jmax);
    if (velmax > 0.5)
      {
        printf("Warning> t=%d velmax=%g x=%d y=%d\n", t, 
               velmax, imax, jmax);
        printf("CFL number too large. Try smaller tau\n");
        exit(1);
      }
  }      
#endif

  /* ---------------------------------------------------------------
     Outputs
     ---------------------------------------------------------------- */
  if (t % out.txplot == 0)
    printxpl(t / out.txplot, u, v, cc, p, 
             num.nx, num.ny, num.tau, num.du, num.dp);
  
  if (t % out.tprint == 0)
    {
      gnuvectors(u, a1, v, a2, num.nx, num.ny, "vectors");

      if (t == 0)
        msgptr = fopen("msg", "wt");
      else
        msgptr = fopen("msg", "at");

      printf("t=%d velmax=%g\n", t, velmax);
      printf("t=%d ncycle=%d npre=%d npost=%d div max before=%g %d %d div mean before=%g\n",
             t , num.cycles, num.npre, num.npost,
             fmodmax(divs, num.nx, num.ny, &imax, &jmax)/h, imax, jmax, 
             fmean(divs, 0.0, 1, num.nx, num.ny)/h);
      printf("t=%d div max=%g %d %d div mean=%g\n", 
             t ,
             fmodmax(work, num.nx, num.ny, &imax, &jmax)/h, imax, jmax, 
             fmean(work, 0.0, 1, num.nx, num.ny)/h);
      fprintf(msgptr, "%d %g %g %g %g %g %g %g %g %g %g\n",
              t,
              phys.t,
              velmax * h / num.tau * num.du,
              interfaces_axivolume(interfaces, Nint),
              findmax(cc, num.nx, num.ny, &imax, &jmax),
              fmodmax(work, num.nx, num.ny, &imax, &jmax)/h,
              num.tau,
              interfaces_x_amp(interfaces, Nint)/(num.nx - 2),
              interfaces_y_amp(interfaces, Nint)/(num.ny - 2),
              pg*num.dp, 
              pe*num.dp);
      printf("mass=%g\n", sumfield(cc, num.nx, num.ny));
      printf("min c: %g  max c: %g\n", findmin(cc, num.nx, num.ny, 
                                               &imax, &jmax),
             findmax(cc, num.nx, num.ny,
                     &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(interfaces[i], msgptr, 100);
      fclose(msgptr);

      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 < Nint; i++) {
        fprintf(msgptr, "\n");
        for (j = 0; j < interfaces[i].n - 1; j++)
          fprintf(msgptr, "%g\n", interfaces[i].t[j]);
      }
      fclose(msgptr);
    }
  
  /* ---------------------------------------------------------------
     Iteration
     ---------------------------------------------------------------- */
interu = fopen("interu", "wt");
interv = fopen("interv", "wt");     
  /*------------------- Advection ------------------- */
  for (i = 0; i < Nint; i++) {
    interface_advect(interfaces[i], u, v, a1, a2, cc, num.nx, num.ny);
    interface_splines(interfaces[i], num.nx, num.ny, AXIS);
#ifndef NOREDIS
    interface_redis(&interfaces[i], num.lredis);
    interface_splines(interfaces[i], num.nx, num.ny, AXIS);
#endif
  }
fclose(interu); fclose(interv);

  interface_to_c(interfaces, Nint,
                 intersection, S12, cc, 0.0, 0.0, num.nx, num.ny,
                 NOTENSION, NULL, NULL, 0.0, 0.0);
  interface_to_c(interfaces, Nint,
                 intersection, S12, a1, 0.5, 0.0, num.nx, num.ny,
                 NOTENSION, NULL, NULL, 0.0, 0.0);
  interface_to_c(interfaces, Nint,
                 intersection, S12, a2, 0.0, 0.5, num.nx, num.ny,
                 NOTENSION, NULL, NULL, 0.0, 0.0);
  bc_scalar(a1, num.nx, num.ny, NULGRAD);
  bc_scalar(a2, num.nx, num.ny, NULGRAD2);
  bc_scalar(cc, num.nx, num.ny, NULGRAD);
}

---