---

ini_conditions.c

Go to the documentation of this file.

/*------------------------------------------------------------------*/
/* SCCS Information: %W%    %G% */
/*------------------------------------------------------------------*/
#include "utilf.h"
#include "stypes.h"
#include <readSIunits.h>
#include "markers.h"
#include "ini_conditions.h"
#include "plgndr.h"

#define YMIN 1.9999


static real angle(real x, real y)
{
  real theta;
  
  if (x == 0.0) {
    if (y >= 0.0)
      return PI/2.;
    return -PI/2.;
  }
  theta = atan(y/x);
  if (x < 0.0)
    theta += PI;
  return theta;
}


static void inidrop(interface *in,
                    real2D u, real2D v, real2D p,
                    real sigma,
                    real tau, 
                    real xvel,
                    real xc, real diameter,
                    real eps, int mode,
                    int nx, int ny)
{
  real h, r, tpi, alpha;
  int i, j;
  
  h = 1.0 / (nx - 2);
  
  printf("------- Initialisation parameters ------\n");
  printf("Droplet: C=(%g,0), D=%g, mode=%d, eps=%g\n", 
         xc, diameter, mode, eps);
    
  diameter *= 0.5 * (nx - 2);
  printf("radius: %f\n", diameter);
  sigma *= tau * tau / (h * h * h);
  xc = xc*(nx - 2) + 2.0;

  for (i = 1; i <= nx; i++)
    for (j = 1; j <= ny; j++) { 
      /*
        u[i][j] = xvel*tau/h*(i - 0.5*(nx + 1))/(0.5*(nx + 1));
        v[i][j] = xvel*tau/h*(j - 0.5*(ny + 1))/(0.5*(nx + 1));
        */
      r = sqrt(sq(i + 0.5 - xc) + sq(j - 1.5));
      alpha = angle(i + 0.5 - xc, j - 1.5);
      if (r <= diameter*(1. + eps*sphericalY(mode, 0, alpha)))
        u[i][j] = xvel*tau/h;
      p[i][j] = - 2.*sigma/diameter;
    }

  in->n = 10*nx;
  in->x = (real *)malloc(in->n * sizeof(real));
  in->y = (real *)malloc(in->n * sizeof(real));
  in->t = (real *)malloc(in->n * sizeof(real));
  in->spx = (polynom3 *)malloc((in->n - 1) * sizeof(polynom3));
  in->spy = (polynom3 *)malloc((in->n - 1) * sizeof(polynom3));
  in->p = (real *)malloc(in->n * sizeof(real));
  
  for (i = 0; i < in->n; i++) {
    tpi = PI*i/(in->n - 1);
    r = diameter*(1. + eps*sphericalY(mode, 0, tpi));
    in->x[i] = xc + r*cos(tpi);
    in->y[i] = YMIN + r*sin(tpi);
  }
  in->bc = AXITRUE;
}


void inidroplet(interface **in,
                int *nin,
                real2D u, real2D v, real2D p,
                real sigma,
                real tau, 
                real xvel,
                int nx, int ny)
{
  real diameter, xc, eps = 0.0;
  int mode = 1;
  
  diameter = 0.5; xc = 0.5;

  realreadSI("diameter", &diameter, stdin);
  realreadSI("xc", &xc, stdin);
  intread("mode", &mode, stdin);
  realreadSI("eps", &eps, stdin);

  *in = (interface *)malloc(sizeof(interface));
  *nin = 1;

  inidrop(*in, u, v, p, sigma, tau, xvel, xc, diameter, eps, mode, nx, ny);
}


void inidroplets(interface **in,
                int *nin,
                real2D u, real2D v, real2D p, 
                real sigma,
                real tau, 
                real xvel,
                int nx, int ny)
{
  real diameter, xc, yc, eps = 0.0;
  int mode = 1;
  
  diameter = 0.5; xc = 0.5; yc = 0.5;

  realreadSI("diameter", &diameter, stdin);
  realreadSI("xc", &xc, stdin);
  realreadSI("yc", &yc, stdin);
  intread("mode", &mode, stdin);
  realreadSI("eps", &eps, stdin);

  *in = (interface *)malloc(2 * sizeof(interface));
  *nin = 2;
  inidrop(&(*in)[0], u, v, p, sigma, tau, xvel, 
          xc - 0.6*diameter, 
          diameter, eps, mode, nx, ny);
  inidrop(&(*in)[1], u, v, p, sigma, tau, -xvel, xc + 0.6*diameter, 
          diameter, eps, mode, nx, ny);
}


void inifusion(interface *in, int nx, int ny)
{
  real h, r, tpi;
  real diameter = 0.4, eps = 0.0, xc1, xc2;
  int mode = 1;
  int i;
  real width = 3.0, height = 2.0, theta;
  
  h = 1.0 / (nx - 2);

  realreadSI("diameter", &diameter, stdin);
  intread("mode", &mode, stdin);
  realreadSI("eps", &eps, stdin);
  
  printf("------- Initialisation parameters ------\n");
  printf("Fusion: D=%g, mode=%d, eps=%g\n", diameter, mode, eps);
    
  diameter *= 0.5 * (nx - 2);
  printf("radius: %f\n", diameter);
  xc1 = 0.5*(nx - 2) + 2.0 + diameter + width/2;
  xc2 = 0.5*(nx - 2) + 2.0 - diameter - width/2;
  theta = (height + width/2.)/diameter;

  in->n = nx*5 + 1;
  in->x = (real *)malloc(in->n * sizeof(real));
  in->y = (real *)malloc(in->n * sizeof(real));
  in->t = (real *)malloc(in->n * sizeof(real));
  in->spx = (polynom3 *)malloc((in->n - 1) * sizeof(polynom3));
  in->spy = (polynom3 *)malloc((in->n - 1) * sizeof(polynom3));
  in->p = (real *)malloc(in->n * sizeof(real));
  
  for (i = 0; i < (in->n - 1)/2; i++) {
    tpi = (PI - theta)*i/((in->n - 1)/2 - 1);
    r = diameter*(1. + eps*sphericalY(mode, 0, tpi));
    in->x[i] = xc1 + r*cos(tpi);
    in->y[i] = YMIN + r*sin(tpi);
  }
  in->x[(in->n - 1)/2] = nx/2 + 1;
  in->y[(in->n - 1)/2] = YMIN + height;
  for (i = (in->n - 1)/2 + 1; i < in->n; i++) {
    tpi = theta + (PI - theta)*(i - (in->n - 1)/2 - 1)/((in->n - 3)/2);
    r = diameter*(1. + eps*sphericalY(mode, 0, tpi));
    in->x[i] = xc2 + r*cos(tpi);
    in->y[i] = YMIN + r*sin(tpi);
  }
  in->bc = AXITRUE;
}


void ini_impact(interface *in, real2D u, real2D v, real xvel, real tau,
                int nx, int ny)
{
  real h, tpi, r;
  real diameter = 0.2, thickness = 0.3, xc1, xc2;
  int i, j;
  real width = 2.0, height = 4.0, theta;
  
  h = 1.0 / (nx - 2);

  realreadSI("diameter", &diameter, stdin);
  realreadSI("thickness", &thickness, stdin);
  
  printf("------- Initialisation parameters ------\n");
  printf("Impact: D=%g thickness=%g xvel=%g\n", diameter, thickness, xvel);
    
  diameter *= 0.5 * (nx - 2);
  printf("radius: %f\n", diameter);
  xc1 = thickness*(nx - 2) + 2.0 + diameter + width/2;
  xc2 = thickness*(nx - 2) + 2.0 - diameter - width/2;
  theta = (height + width/2.)/diameter;

#if 0
  for (i = 2; i < nx; i++)
    for (j = 2; j < ny; j++) {
      r = sqrt(sq((real)i - xc1) + sq((real)j - 1.5));
      if (r <= diameter)
        u[i][j] = xvel*tau/h;
    }
#else
  for (i = 2 /* xc1 - diameter */ ; i < nx; i++)
    for (j = 2; j < ny; j++)
        u[i][j] = xvel*tau/h;
#endif

  in->n = nx*5 + 1;
  in->x = (real *)malloc(in->n * sizeof(real));
  in->y = (real *)malloc(in->n * sizeof(real));
  in->t = (real *)malloc(in->n * sizeof(real));
  in->spx = (polynom3 *)malloc((in->n - 1) * sizeof(polynom3));
  in->spy = (polynom3 *)malloc((in->n - 1) * sizeof(polynom3));
  in->p = (real *)malloc(in->n * sizeof(real));
  
  for (i = 0; i < (in->n - 1)/2; i++) {
    tpi = (PI - theta)*i/((in->n - 1)/2 - 1);
    in->x[i] = xc1 + diameter*cos(tpi);
    in->y[i] = YMIN + diameter*sin(tpi);
  }
  in->x[(in->n - 1)/2] = (xc1 + xc2)/2.;
  in->y[(in->n - 1)/2] = YMIN + height;
  for (i = (in->n - 1)/2 + 1; i < in->n; i++) {
    in->x[i] = xc2 + diameter;
    in->y[i] = YMIN + height + width/2. + (i - (in->n - 1)/2 - 1)*((real)ny - 2. - height - width/2.)/(real)(in->n - 2 - (in->n - 1)/2);
  }
  in->y[in->n - 1] = 0.001 + (real)ny;
  in->bc = AXIBOTH;
#if 0
  for (i = 0; i < in->n/2; i++) {
    double tmp;
    tmp = in->x[i];
    in->x[i] = in->x[in->n - 1 - i];
    in->x[in->n - 1 - i] = tmp;

    tmp = in->y[i];
    in->y[i] = in->y[in->n - 1 - i];
    in->y[in->n - 1 - i] = tmp;
  }
#endif
#if 0  
  for (i=0; i<=in->n-1; i++)
          in->x[i] = (real)nx - in->x[i] + 0.05;
#endif
}


void inifile(interface *in, int nx, int ny, char *file)
{
  FILE *fptr = fopen(file, "rt");
  int i, n = 0;
  float x, y;

  if (!fptr) {
    fprintf(stderr, "inifile(): cannot open file `%s'\n", file);
    exit(1);
  }
  while (fscanf(fptr, "%f %f", &x, &y) == 2)
    n++;
  printf("------- Initialisation parameters ------\n");
  printf("File: %s, npoints: %d\n", file, n);
  in->n = n;
  in->x = (real *)malloc(in->n * sizeof(real));
  in->y = (real *)malloc(in->n * sizeof(real));
  in->t = (real *)malloc(in->n * sizeof(real));
  in->spx = (polynom3 *)malloc((in->n - 1) * sizeof(polynom3));
  in->spy = (polynom3 *)malloc((in->n - 1) * sizeof(polynom3));
  in->p = (real *)malloc(in->n * sizeof(real));
  
  rewind(fptr);
  for (i = 0; i < n; i++) {
    fscanf(fptr, "%f %f", &x, &y);
    in->x[i] = (real)(nx - 2)*x + 2.;
    in->y[i] = (real)(nx - 2)*y + YMIN;
  }
/*cld*/  
  in->bc = AXIBOTH;
/*cld*/
/*in->bc = AXITRUE*/
  fclose(fptr);
}


static void inibub(interface *in,
                   real2D u, real2D v, real2D p, 
                   real sigma,
                   real tau, 
                   real xvel,
                   real po,
                   real xc, real diameter,
                   real eps, int mode,
                   int nx, int ny)
{
  real h, r, tpi;
  int i, j;
  
  h = 1.0 / (nx - 2);
  
  printf("------- Initialisation parameters ------\n");
  printf("Bubble: C=(%g,0), d=%g, mode=%d, eps=%g\n", 
         xc, diameter, mode, eps);
    
  diameter *= 0.5 * (nx - 2);
  printf("radius: %f\n", diameter);
  sigma *= tau * tau / (h * h * h);
  po *= sq(tau/h);
  xc = xc*(nx - 2) + 2.0;

  for (i = 1; i <= nx; i++)
    for (j = 1; j <= ny; j++) { 
      /*
        u[i][j] = xvel*tau/h*(i - 0.5*(nx + 1))/(0.5*(nx + 1));
        v[i][j] = xvel*tau/h*(j - 0.5*(ny + 1))/(0.5*(nx + 1));

        r = sqrt(sq(i - xc) + sq(j - 1.5));
        u[i][j] = -sq(diameter/r)/r*((real)i - xc)*xvel*tau/h;
        r = sqrt(sq(i + 0.5 - xc) + sq(j - 2.0));
        v[i][j] = -sq(diameter/r)/r*((real)j - 2.0)*xvel*tau/h;
      */
      p[i][j] = - 2.*sigma/diameter; 
    }
  
  in->n = 10*nx;
  in->x = (real *)malloc(in->n * sizeof(real));
  in->y = (real *)malloc(in->n * sizeof(real));
  in->t = (real *)malloc(in->n * sizeof(real));
  in->spx = (polynom3 *)malloc((in->n - 1) * sizeof(polynom3));
  in->spy = (polynom3 *)malloc((in->n - 1) * sizeof(polynom3));
  in->p = (real *)malloc(in->n * sizeof(real));
  
  for (i = 0; i < in->n; i++) {
    tpi = PI*(in->n - 1 - i)/(in->n - 1);
    r = diameter*(1. + eps*sphericalY(mode, 0, tpi));
    in->x[i] = xc + r*cos(tpi);
    in->y[i] = YMIN + r*sin(tpi);
  }
  in->bc = AXITRUE;
}


void inibubble(interface **in,
               int *nin,
               real2D u, real2D v, real2D p, 
               real sigma,
               real tau, 
               real xvel,
               real po,
               int nx, int ny)
{
  real diameter, xc, eps = 0.0;
  int mode = 1;
  
  diameter = 0.5; xc = 0.5;

  realreadSI("diameter", &diameter, stdin);
  realreadSI("xc", &xc, stdin);
  intread("mode", &mode, stdin);
  realreadSI("eps", &eps, stdin);

  *in = (interface *)malloc(sizeof(interface));
  *nin = 1;

  inibub(*in, u, v, p, sigma, tau, xvel, po, xc, diameter, eps, mode, nx, ny);
}


void iniwave(interface **in,
             int *nin,
             real2D u, real2D v, real2D p, 
             real tau, 
             real xvel,
             int nx, int ny)
{
  real h, x, y, ampli = 0.1, wavenr = 1.0, position = 0.5;
  int i;

  realreadSI("amplitude", &ampli, stdin);
  realreadSI("wavenr", &wavenr, stdin);
  realreadSI("position", &position, stdin);
  printf("------- Initialisation parameters ------\n");
  printf("Wave: wavenr=%g amplitude=%g position=%g\n", 
         wavenr, ampli, position);
  
  ampli *= ny - 2;
  h = 1.0 / (nx - 2);

  *in = (interface *)malloc(sizeof(interface));
  *nin = 1;
  (*in)->n = 10*nx;
  (*in)->x = (real *)malloc((*in)->n * sizeof(real));
  (*in)->y = (real *)malloc((*in)->n * sizeof(real));
  (*in)->t = (real *)malloc((*in)->n * sizeof(real));
  (*in)->spx = (polynom3 *)malloc(((*in)->n - 1) * sizeof(polynom3));
  (*in)->spy = (polynom3 *)malloc(((*in)->n - 1) * sizeof(polynom3));
  (*in)->p = (real *)malloc((*in)->n * sizeof(real));
  
  for (i = 0; i < (*in)->n - 1; i++) {
    x = (real)i * (real)(nx - 2.0) / (real)((*in)->n - 1) + 2.1;
    y = ampli*cos((x - 2.0)*2.0*PI*wavenr/(real)(nx - 2.0)) 
      + position*ny + 1.501;
    (*in)->x[i] = x;
    (*in)->y[i] = y;
  }  
  (*in)->x[(*in)->n - 1] = (*in)->x[0] + nx - 2.0;
  (*in)->y[(*in)->n - 1] = (*in)->y[0];
}


void inirealwave(interface **in,
                 int *nin,
                 real2D u, real2D v, real2D p, 
                 real tau, 
                 real xvel,
                 int nx, int ny)
{
  real h, x, y, ampli = 0.1, wavenr = 1.0;
  int i;

  realreadSI("amplitude", &ampli, stdin);
  realreadSI("wavenr", &wavenr, stdin);
  printf("------- Initialisation parameters ------\n");
  printf("Wave: wavenr=%g amplitude=%g\n", wavenr, ampli);
  
  ampli *= ny - 2;
  h = 1.0 / (nx - 2);

  *in = (interface *)malloc(sizeof(interface));
  *nin = 1;
  (*in)->n = 100;
  (*in)->x = (real *)malloc((*in)->n * sizeof(real));
  (*in)->y = (real *)malloc((*in)->n * sizeof(real));
  (*in)->t = (real *)malloc((*in)->n * sizeof(real));
  (*in)->spx = (polynom3 *)malloc(((*in)->n - 1) * sizeof(polynom3));
  (*in)->spy = (polynom3 *)malloc(((*in)->n - 1) * sizeof(polynom3));
  (*in)->p = (real *)malloc((*in)->n * sizeof(real));
  
  for (i = 0; i < (*in)->n - 1; i++) {
    x = i * (nx - 2.0) / (real)((*in)->n - 1) + 2.1;
    y = ampli*cos((x - 2.0)*2.0*PI*wavenr/(nx - 2.0)) + 0.5*ny + 1.001;
    (*in)->x[i] = x;
    (*in)->y[i] = y;
  }  
  (*in)->x[(*in)->n - 1] = (*in)->x[0] + nx - 2.0;
  (*in)->y[(*in)->n - 1] = (*in)->y[0];
}


void initorus(interface *in,
              real2D u, real2D v, real2D p, 
              real tau, 
              real xvel,
              int nx, int ny)
{
  real h, r, tpi;
  real xc, yc, diameter, eps;
  int mode;
  int i, j;
  
  h = 1.0 / (nx - 2);

  realreadSI("diameter", &diameter, stdin);
  realreadSI("xc", &xc, stdin);
  realreadSI("yc", &yc, stdin);
  intread("mode", &mode, stdin);
  realreadSI("eps", &eps, stdin);

  printf("------- Initialisation parameters ------\n");
  printf("torus: C=(%g,%g), d=%g, mode=%d, eps=%g\n", 
         xc, yc, diameter, mode, eps);
  
  diameter *= 0.5 * (nx - 2);
  printf("radius: %f\n", diameter);
  xc = xc*(nx - 2) + 2.0;
  yc = yc*(ny - 2) + 2.0;

  in->n = 10*nx;
  in->x = (real *)malloc(in->n * sizeof(real));
  in->y = (real *)malloc(in->n * sizeof(real));
  in->t = (real *)malloc(in->n * sizeof(real));
  in->spx = (polynom3 *)malloc((in->n - 1) * sizeof(polynom3));
  in->spy = (polynom3 *)malloc((in->n - 1) * sizeof(polynom3));
  in->p = (real *)malloc(in->n * sizeof(real));

  for (i = 0; i < in->n; i++) {
    tpi = 2.*PI*(in->n - 1 - i)/(in->n - 1);
    r = diameter*(1. + eps*sphericalY(mode, 0, tpi));
    in->x[i] = xc + r*cos(tpi);
    in->y[i] = yc + r*sin(tpi);
  }
  in->x[in->n - 1] = in->x[0];
  in->y[in->n - 1] = in->y[0];
  in->bc = TRUE_PERIODIC;
}


void inijet(interface *in,
            real2D u, real2D v, real2D p, 
            real tau, 
            real xvel,
            int nx, int ny)
{
  real h = 1.0/(nx - 2);
  int i, j;

  for (i = nx/3; i <= nx; i++)
    for (j = 2; j < 10; j++)
      u[i][j] = xvel*tau/h*sq(10 - j)/sq(8);
  
  in->n = 4;
  in->x = (real *)malloc(in->n * sizeof(real));
  in->y = (real *)malloc(in->n * sizeof(real));
  in->t = (real *)malloc(in->n * sizeof(real));
  in->spx = (polynom3 *)malloc((in->n - 1) * sizeof(polynom3));
  in->spy = (polynom3 *)malloc((in->n - 1) * sizeof(polynom3));
  in->p = (real *)malloc(in->n * sizeof(real));

  in->x[0] = (nx - 2)/4. + 2.2; in->y[0] = YMIN;
  in->x[1] = (nx - 2)/2. + 2.; in->y[1] = (ny - 2)/2. + 2.;
  in->x[2] = 3.*(nx - 2)/4. + 2.; in->y[2] = (ny - 2)/4. + 2.;
  in->x[3] = (nx - 2)/3. + 2.; in->y[3] = YMIN;
  in->bc = AXITRUE;
}











void ini_impact2(interface *in, real2D u, real2D v, real xvel, real tau,
                int nx, int ny)
{
  real h, tpi, r;
  real diameter = 0.2, thickness = 0.3, xc1, xc2;
  int i, j;
  real width = 4.0, height = 4.0, theta;
  
  h = 1.0 / (nx - 2);

  realreadSI("diameter", &diameter, stdin);
  realreadSI("thickness", &thickness, stdin);
  
  printf("------- Initialisation parameters ------\n");
  printf("Impact: D=%g thickness=%g xvel=%g\n", diameter, thickness, xvel);
    
  diameter *= 0.5 * (nx - 2);
  printf("radius: %f\n", diameter);
  xc1 = thickness*(nx - 2) + 2.0 + diameter + width/2;
  xc2 = thickness*(nx - 2) + 2.0 - diameter - width/2;
  theta = (height + width/2.)/diameter;

#if 0
  for (i = 2; i < nx; i++)
    for (j = 2; j < ny; j++) {
      r = sqrt(sq((real)i - xc1) + sq((real)j - 1.5));
      if (r <= diameter)
        u[i][j] = xvel*tau/h;
    }
#else
  for (i = 2 /* xc1 - diameter */ ; i < nx; i++)
    for (j = 2; j < ny; j++)
        u[i][j] = xvel*tau/h;
#endif

  in->n = nx*5 + 1;
  in->x = (real *)malloc(in->n * sizeof(real));
  in->y = (real *)malloc(in->n * sizeof(real));
  in->t = (real *)malloc(in->n * sizeof(real));
  in->spx = (polynom3 *)malloc((in->n - 1) * sizeof(polynom3));
  in->spy = (polynom3 *)malloc((in->n - 1) * sizeof(polynom3));
  in->p = (real *)malloc(in->n * sizeof(real));
  
  for (i = 0; i < (in->n - 1)/2; i++) {
    tpi = (PI - theta)*i/((in->n - 1)/2 - 1);
    in->x[i] = xc1 + diameter*cos(tpi);
    in->y[i] = YMIN + diameter*sin(tpi);
  }
  in->x[(in->n - 1)/2] = (xc1 + xc2)/2.;
  in->y[(in->n - 1)/2] = YMIN + height;
  for (i = (in->n - 1)/2 + 1; i < in->n; i++) {
    in->x[i] = xc2 + diameter;
    in->y[i] = YMIN + height + width/2. + (i - (in->n - 1)/2 - 1)*((real)ny - 2. - height - width/2.)/(real)(in->n - 2 - (in->n - 1)/2);
  }
  in->y[in->n - 1] = 0.001 + (real)ny;
  in->bc = AXIBOTH;
#if 0
  for (i = 0; i < in->n/2; i++) {
    double tmp;
    tmp = in->x[i];
    in->x[i] = in->x[in->n - 1 - i];
    in->x[in->n - 1 - i] = tmp;

    tmp = in->y[i];
    in->y[i] = in->y[in->n - 1 - i];
    in->y[in->n - 1 - i] = tmp;
  }
#endif
#if 1  
  for (i=0; i<=in->n-1; i++)
          in->x[i] = (real)nx - in->x[i] + 0.05;
#endif
}

---