---

Shear.c

Go to the documentation of this file.

/*-----------------------------------------------------------------*/
/* SCCS Information: %W%    %G% */
/*-----------------------------------------------------------------*/
#include "utilf.h"
#include "markers.h"
#include "surfaces.h"
#include "extra.h"
#include "global.h"
#include "ini_conditions.h"
#include "save.h"
#include "make_bc.h"
#include "mg.h"
#include "save.h"
#include <readSIunits.h>
#include <time.h>


void initialize()
{
  int restart, i, j;
  FILE *fptr, *fptr1, *fptr2;

  readglob();
  initglobal();

  restart = 0;
  intread("restart", &restart, stdin);
  printf("restart=%d\n", restart);
  
  printhist();
  printnumbers();

  switch(num.ng = computeng(num.nx, num.ny)) {
  case NXERROR:
    printf("computeng(): bad nx size\n");
    exit(NXERROR);
  case NYERROR:
    printf("computeng(): bad ny size\n");
    exit(NYERROR);
  case NGERROR:
    printf("computeng(): bad ng size\n");
    exit(NGERROR);
  }
  if (initpressure(num.nx, num.ny, num.ng)) {
    printf("initpressure(): cannot allocate memory\n");
    exit(1);
  }

  if (restart == 0) {
    interfaces = (interface *)malloc(sizeof(interface));
    Nint = 1;
        inifile(interfaces, num.nx, num.ny, "file.init"); 
      //  inifusion(interfaces, num.nx, num.ny); 
//    ini_impact(interfaces, u, v, phys.xvel, num.tau, num.nx, num.ny); 
//      ini_impact2(interfaces, u, v, phys.xvel, num.tau, num.nx, num.ny);
/*        inidroplet(&interfaces, &Nint, 
          u, v, p, 
          phys.sigma, num.tau, phys.xvel,
          num.nx, num.ny); */
    /* inibubble(&interfaces, &Nint, 
              u, v, p, 
              phys.sigma, num.tau, phys.xvel, phys.po,
              num.nx, num.ny); */
    /*
      iniwave(&interfaces, &Nint, 
      u, v, p, num.tau, phys.xvel, num.nx, num.ny);
      */
    /*
      interfaces = (interface *)malloc(sizeof(interface));
      Nint = 1;
      inijet(interfaces, u, v, p, num.tau, phys.xvel, num.nx, num.ny);
    */
  }
  else if (restart == 1) {
    iniget(u, v, ap, p, &interfaces, &Nint, &phys.po, &phys.vo, 
           &phys.t, &phys.pext, 
           &num.tau, num.nx, num.ny, "save.dat");
    printf("fields u v ap p read\n");
  }

  fptr = fopen("markers.init", "wt");
  fptr1 = fopen("curvature.init", "wt");
  fptr2 = fopen("splines.init", "wt");
  for (i = 0; i < Nint; i++) {
    interface_splines(interfaces[i], num.nx, num.ny);
    interface_write_markers(interfaces[i], fptr);
    interface_write_splines(interfaces[i], fptr2, 1000);
    curvature_write(interfaces[i], fptr1);
    interface_redis(&interfaces[i], num.lredis/2.);
    interface_splines(interfaces[i], num.nx, num.ny);
  }
  fclose(fptr); fclose(fptr1); fclose(fptr2);

  printf("Filling fractions ... ");
  interface_fill(interfaces, 1, ap, num.nx, num.ny);
  printf("Ok\n");
  
  printf("Initializing fraction fields ... "); fflush(stdout);
  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);

  for (i = 2; i <= num.nx; i++)
    for (j = 2; j <= num.ny; j++) {
      au[i][j] = (ap[i][j] + ap[i-1][j])/2.;
      av[i][j] = (ap[i][j] + ap[i][j-1])/2.;
    }
  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);
  printf("Ok\n");

  /*
    if (restart == 0)
    phys.vo = fabs(interfaces_axivolume(interfaces, Nint));
  */  
  bc_vector_bound(u, v, num.nx, num.ny);

  printxpl(u, v, ap, p, num.nx, num.ny, num.tau, num.du, num.dp, "init.datb");

}


int main(int argc, char *argv[])
{
  int t, steptime;
  long CPUtime;
  
  initialize();
  CPUtime = 0;
  for (t = 0; t <= num.tmax; t++) {
    steptime = clock();
    if (t > 0 && t % out.tsave == 0) {
      rename("save.dat", "save.dat.bak");
      save(u, v, ap, p, interfaces, Nint, phys.po, phys.vo, phys.t, phys.pext, 
           num.tau, num.nx, num.ny, "save.dat");
      printf("fields u v ap p written\n");
    }
    timestep(t);
    CPUtime += clock() - steptime ;
  }
  printf("CPU time: %g s\n", (float) CPUtime / 
         (float) CLOCKS_PER_SEC);
  printf("Grid points / sec: %g\n",(float) (num.tmax + 1) * num.nx * num.ny
         / (float) CPUtime * CLOCKS_PER_SEC);
}

---