sandbox/yonghui/archrrc.c
real 3d arch with KV model ADD THE RRC LATER
#include "grid/octree.h"
#include "navier-stokes/centered.h"
#include "utils.h"
#include "distance.h"
#include "fractions.h"
#include "view.h"
#include "lambda2.h"
#define MIN_LEVEL 5
#define LEVEL 7
#define MAX_LEVEL 10
#define eq_r 1.
#define ADAPT 0
#define tmax 1*2.*M_PI
#define alpha_w 10.
Importing the geometry
Importing the geometry
void fraction_from_stl (scalar f, FILE * fp)
{
We read the STL file and compute the bounding box of the model.
coord * p = input_stl (fp);
coord min, max;
bounding_box (p, &min, &max);
double maxl = -HUGE;
foreach_dimension()
if (max.x - min.x > maxl)
maxl = max.x - min.x;
scalar d[];
distance (d, p);
vertex scalar phi[];
foreach_vertex()
phi[] = (d[] + d[-1] + d[0,-1] + d[-1,-1] +
d[0,0,-1] + d[-1,0,-1] + d[0,-1,-1] + d[-1,-1,-1])/8.;
fractions (phi, f);
}
volume fraction field declaire
scalar f0[];
int main()
{
We read the STL file, compute the bounding box of the model and set the domain center and size using this bounding box.
init_grid (8);
size (8.);
origin (-1*L0/4 ,-2*L0/5, 0.);
run();
}
We then set an inflow velocity on the top and free outflow on the bottum. axi- Y Z - front z =L0 back z=0 Y - top bottom X - left right
u.n[back] = ( y >= 0)? dirichlet(10.*f0[]*cos(t)):neumann(0);
u.t[back] = ( y >= 0)? dirichlet(0):neumann(0);
u.r[back] = ( y >= 0)? dirichlet(0):neumann(0);
p[back] = ( y >= 0)? dirichlet(0):neumann(0) ;
pf[back] = ( y >= 0)? dirichlet(0):neumann(0) ;
p[front] = neumann(0) ;
pf[front] = neumann(0) ;
event init (t = 0) {
if (!restore (file = "restart")) {
FILE * fp = fopen ("aortic_arch.stl", "r");
char name2[10];
for (int ii=MIN_LEVEL; ii<=LEVEL; ii++){
fraction_from_stl (f0, fp);
adapt_wavelet ({f0}, (double[]){0.00001}, ii);
//dump file
sprintf (name2, "dump_%d",ii);
dump(file = name2);
}
fraction_from_stl (f0, fp);
view (fov = 19.9677, quat = {0.614306,0.390179,0.155159,0.66807},
tx = -0.120804, ty = -0.340923, width = 1920, height = 1080);
draw_vof ("f0");
squares ("f0", min = -10, max =10);
cells();
save ("init.png");
double viscosity = sq(eq_r) / sq(alpha_w);
const face vector muc[] = {viscosity, viscosity, viscosity};
mu = muc;
}
}
event bc (t <= tmax; i += 1) {
//BC : lateral velocity = 0
foreach()
foreach_dimension()
u.x[] = f0[]*u.x[];
boundary ((scalar *){u , p});
//position du test point
double px = 3.5;
double py = 0.8;
double pz = 2.;
fprintf(stderr, "%d %g %g %g \n", i, t, dt,interpolate(u.x , px, py, pz))
}
event snapshot (i += 1000)
{
char name[80];
sprintf (name, "dump-%d", i);
scalar l2[];
lambda2 (u,l2);
dump (file = name);
}
int j = 0;
event profill (t += 0.4*M_PI)
{
if (t > tmax - 2.*M_PI){
j += 1;
char name2[80];
sprintf (name2, "test-%d", j);
FILE *fp2;
fp2 = fopen (name2, "w");
double uyy = -0.5 ;
double uxx = 2.5;
double velou,velov,velo;
for (double uzz = 0.5; uzz <= 4.; uzz += 4./100. ){
velou = interpolate(u.x , uxx, uyy, uzz);
velov = interpolate(u.y , uxx, uyy, uzz);
velo = velou * cos(M_PI/4.) + velov * sin(M_PI/4.);
fprintf(fp2,"%g %g %g %g %g\n", t, uzz-2.5, velo, velou, velov);
}
}
}
event video (t += tmax/300.){
scalar l2[];
lambda2 (u,l2);
input and output view
view (fov = 20, quat = {0,0,0,1}, tx = -0.19684, ty = 0.00981146, bg = {0.3,0.4,0.6}, width = 640, height = 640);
// draw_vof("f0");
squares ("u.z", alpha = 0.4);
cells( alpha = 0.4);
save ("input_uz.mp4");
main chanel view
clear();
view (fov = 18.0787, quat = {0.126228,-0.561414,-0.620277,0.533052}, tx = 0.00499754, ty = -0.373772, bg = {0.3,0.4,0.6}, width = 640, height = 640, samples = 1);
//draw_vof("f0");
squares("u.z",n = {1,-1,0},alpha =2.4 );
squares("u.z",n = {-1,-1,1.5},alpha = -3.5 );
cells(n = {1,-1,0},alpha =2.4 );
cells(n = {-1,-1,1.5},alpha = -3.5 );
save ("mainchannel_uz.mp4");
clear();
view (fov = 18.0787, quat = {0.126228,-0.561414,-0.620277,0.533052}, tx = 0.00499754, ty = -0.373772, bg = {0.3,0.4,0.6}, width = 640, height = 640, samples = 1);
//draw_vof("f0");
squares("l2",n = {1,-1,0},alpha =2.4 );
squares("l2",n = {-1,-1,1.5},alpha = -3.5 );
cells(n = {1,-1,0},alpha =2.4 );
cells(n = {-1,-1,1.5},alpha = -3.5 );
save ("mainchannel_l2.mp4");
clear();
view (fov = 18.0787, quat = {0.126228,-0.561414,-0.620277,0.533052}, tx = 0.00499754, ty = -0.373772, bg = {0.3,0.4,0.6}, width = 640, height = 640, samples = 1);
//draw_vof("f0");
squares("u.x",n = {1,-1,0},alpha =2.4 );
squares("u.x",n = {-1,-1,1.5},alpha = -3.5 );
cells(n = {1,-1,0},alpha =2.4 );
cells(n = {-1,-1,1.5},alpha = -3.5 );
save ("mainchannel_ux.mp4");
cross section view
view (fov = 11.1239, quat = {0.658889,-0.274942,-0.26525,0.648006}, tx = -0.207846, ty = -0.306235, bg = {0.3,0.4,0.6}, width = 640, height = 640, samples = 1);
// draw_vof("f0");
squares("u.x",n = {-1,-1,0},alpha =-2.4);
cells(n = {-1,-1,0},alpha =-2.4);
save ("section_ux.mp4");
squares("l2",n = {-1,-1,0},alpha =-2.4);
cells(n = {-1,-1,0},alpha =-2.4);
save ("section_l2.mp4");
}
#if ADAPT
event adapt (i++) {
double uemax = 0.01;
adapt_wavelet ({u}, (double[]){uemax,uemax,uemax}, MAX_LEVEL, 4);
}
#endif