sandbox/msaini/tests/embedpotential.c

    I want to share the results from Saini2022 for future references. The main finding is that the potential flow solution for bubble collapse is singular at the contact line. The potential flow solution is given by the solution of 2p=0\displaystyle \boldsymbol{\nabla}^2 p = 0 is domain with BCs as shown in following figure
    Domain and BCs for problem. Definition of different coordinate systems is also shown.

    We use embed boundary to impose boundary conditions the bubble interface and solve the Laplace equation using the Poisson solver.

    #include "embed.h"
    #include "poisson.h"
    #include "utils.h"
    #include "output.h"
    #include "view.h"
    
    double R = 1;
    double pinf = 10.;
    scalar psi[], src[];
    scalar psi1[],psi1b[];
    scalar pr[];
    
    psi[top] = dirichlet (pinf);
    psi[right] = dirichlet (pinf);
    psi[bottom] = neumann(0.);
    psi[left] = neumann(0.);
    
    FILE *fp,*fp2;
    
    int main(int argc, char *argv[]){
    
      int maxlevel = 12;
    
      for(double i = 7.5; i <= 12; i+= 1.5){
        double theta0 = i*pi/18.;
        double xc = cos(theta0);
    
        L0 = 100.;
        size (L0);
        Y0 = 0.;
        X0 = 0.;
      
        init_grid (1 << 8);
        refine ( level <= (10 - (sqrt(x*x + y*y + z*z)/(3.) > 1)*30*pow(sqrt(x*x + y*y + z*z)/(3.) - 1.,1) ));
        refine(level < maxlevel && sqrt(sq(x - xc) + sq(y)+ z*z) > (1. - 10.*sqrt(2.)*100./(1<<12)) && sqrt(sq(x - xc) + sq(y)+ z*z) < (1. + 10.*sqrt(2.)*100./(1<<12)));
      
        vertex scalar phi[];
        foreach_vertex() {
          phi[] = sq(x - xc) + y*y - R*R;
        }
        fractions (phi, cs, fs);
        
    #if TREE
        cs.refine = cs.prolongation = fraction_refine;
    #endif
        restriction ({cs,fs});
        
        cm = cs;
        fm = fs;
        
        foreach (){
          src[] = 0.;
          psi[] = pinf;
        }
          
        psi[embed] = dirichlet(1.);
        
        psi.third = true;
        
        face vector alphav[];
        
        foreach_face (x) {
          alphav.x[] = y*fs.x[];
        }
        
        foreach_face (y) {
          alphav.y[] = y*fs.y[];
        }
        
        poisson (psi, src, alphav, tolerance = 1e-7, minlevel = 4);
    
        scalar gradpn[];
    
        char dpname[50];
        sprintf(dpname,"dpressure-%2.2f.dat",i);
        fp = fopen(dpname,"w");
    
        char dpname2[50];
        sprintf(dpname2,"dpint-%2.2f.dat",i);
        fp2 = fopen(dpname2,"w");
    
        foreach(){
          if (cs[] > 0. && cs[] < 1.){
    	coord n = facet_normal (point, cs, fs), p;
    	double alpha = plane_alpha (cs[], n);
    	double area = plane_area_center (n, alpha, &p);
    	bool dirichlet;
    	double vb = psi.boundary[embed] (point, point, psi, &dirichlet);
    	if (dirichlet){
    	  double val;
    	  normalize (&n);
    	  gradpn[] = dirichlet_gradient (point, psi, cs, n, p, vb, &val)/(1. - pinf);
    	  fprintf(fp,"%10.9f %10.9f %10.9f\n",x,y,fabs(gradpn[]));
    	  fprintf(fp2,"%10.9f %10.9f %10.9f %10.9f\n",atan2(x,y),x,y,fabs(gradpn[]));
    	}
    	else
    	  gradpn[] = 0;
          }
          else{
    	gradpn[] = cs[] >= 1. ? sqrt(sq((psi[1] - psi[-1])/2/Delta) + sq((psi[0,1] - psi[0,-1])/2/Delta))/(1. - pinf): 0.;
    	if(x < 5. && y < 5.)
    	  fprintf(fp,"%10.9f %10.9f %10.9f\n",x,y,fabs(gradpn[]));
          }
        }
    
        view(fov = 4., quat = {0,0,-cos(pi/4.),cos(pi/4.)}, width = 1980, height = 1980);
        box (notics=true);
        draw_vof (c = "cs",lw=4);
        isoline (phi = "gradpn", n = 200, min = -9, max = 10);
        mirror({0,1}) {
          box (notics=true);
          draw_vof("f",lw=4);
          isoline (phi = "gradpn", n = 200, min = -9, max = 10);
        }
        char filename[80];
        sprintf(filename,"isolines-%dpiby18.png",(int)i);
        save(filename);

    We output data in a file to launch a DNS simulation from initial condition given by the potential flow solution. In case of MPI, we write the data in to seperate files that can be combined later by Bash commands for ex. cat.

        char namew[50];
    #if _MPI
        sprintf(namew,"data-%2.2f-%d-%d.dat",theta0,maxlevel,pid());
    #else
        sprintf(namew,"data-%2.2f.dat",theta0);
    #endif
        fp = fopen(namew,"w");
        if(pid() == 0)
          fprintf(fp,"%ld ",grid->tn);
        
        foreach()
          fprintf(fp,"%d %d %d %lf ",point.level,point.i,point.j,psi[]);
        fclose(fp);
        
      }
    
    }
    Shapes at tend (script)

    Shapes at tend (script)

    Contour plots (script)