Commit 751e881f authored by Michel Henry's avatar Michel Henry
Browse files

robust adapt + correct leidenfrost

parent e9c3c9cf
......@@ -1578,7 +1578,10 @@ void fluid_problem_adapt_mesh(FluidProblem *problem, Mesh *new_mesh, int old_n_p
int eid=-1;
double xi_old[D];
mesh_tree_particle_to_mesh(get_mesh_tree(problem), (size_t) 1, x, &eid, xi_old);
if (eid == -1){
printf("No eid found !\n");
exit(0);
}
double s = 0;
double sf_old[fields->element[ifield]->nlocal];
fields->element[ifield]->f(xi_old, sf_old);
......@@ -1657,104 +1660,6 @@ void fluid_problem_adapt_mesh(FluidProblem *problem, Mesh *new_mesh, int old_n_p
}
#if 0
void fluid_problem_adapt_mesh(FluidProblem *problem, Mesh *new_mesh, int old_n_particles, double *old_particle_position, double *old_particle_volume)
{
const FEElement *element = problem->fields->element[0];
int n_fields = fluid_problem_n_fields(problem);
int ndof_new = -1;
for(int c = 0; c < new_mesh->n_elements*problem->field->local_size; ++c){
ndof_new = max(ndof_new, new_map[c]);
}
if(ndof_new < 0){
printf("Invalid parameters !\n");
exit(1);
}
ndof_new++;
double *new_solution = malloc(sizeof(double)*ndof_new);
double *new_xi = malloc(sizeof(double)*new_mesh->n_nodes*D);
int *new_eid = malloc(sizeof(int)*new_mesh->n_nodes);
for(int i = 0; i < new_mesh->n_nodes; ++i)
new_eid = -1;
double *newx = malloc(double)*D*new_mesh->n_nodes);
for (int i = 0;i < new_mesh->n_nodes;++i)
for (int k = 0; k < D; ++k)
newx[i*D+k] = new_mesh->x[i*3+k];
mesh_tree_particle_to_mesh(get_mesh_tree(problem),new_mesh->n_nodes, newx, new_eid, new_xi);
if(problem->n_fluids == 2){
double *concentration_cg = malloc(sizeof(double)*new_mesh->n_elements*N_N);
}
double *new_solution = (double*)malloc(sizeof(double)*new_mesh->n_nodes*n_fields);
double *concentration_cg = (double*)malloc(sizeof(double)*new_mesh->n_elements*N_N);
double *new_xi = (double*)malloc(sizeof(double)*new_mesh->n_nodes*D);
int *new_eid = (int*)malloc(sizeof(int)*new_mesh->n_nodes);
for (int i = 0; i < new_mesh->n_nodes; ++i)
new_eid[i] = -1;
double *newx = (double*)malloc(sizeof(double)*D*new_mesh->n_nodes);
for (int i = 0;i < new_mesh->n_nodes;++i)
for (int k = 0; k < D; ++k)
newx[i*D+k] = new_mesh->x[i*3+k];
mesh_tree_particle_to_mesh(get_mesh_tree(problem),new_mesh->n_nodes, newx, new_eid, new_xi);
for (int i = 0; i < new_mesh->n_nodes; ++i) {
const int *el = problem->mesh->elements+new_eid[i]*N_N;
if(new_eid[i] < 0) {
printf("new mesh point not found in old mesh\n");
exit(1);
}
double phi[element->nlocal];
element->f(new_xi+i*D, phi);
for (int j=0; j<n_fields; ++j) {
new_solution[i*n_fields+j] = 0;
for (int k = 0; k < element->nlocal; ++k)
new_solution[i*n_fields+j] += problem->solution[el[k]*n_fields+j]*phi[k]; // pas juste
}
if (problem->n_fluids == 2) {
concentration_cg[i] = 0;
for (int k = 0; k < element->nlocal; ++k){
concentration_cg[i] += problem->concentration[new_eid[i]*N_N+k]*phi[k];
}
}
}
free(new_eid);
free(new_xi);
free(newx);
fluid_problem_set_mesh(problem,new_mesh);
free(problem->solution);
if (problem->n_fluids == 2) {
free(problem->concentration);
problem->concentration = malloc(sizeof(double)*new_mesh->n_elements*N_N);
for (int iel = 0; iel < new_mesh->n_elements; ++iel) {
const int *el = &new_mesh->elements[iel*N_N];
for (int i = 0; i < N_N; ++i) {
problem->concentration[iel*N_N+i] = concentration_cg[el[i]];
}
}
free(concentration_cg);
}
problem->solution = new_solution;
for (int i = 0; i < problem->n_particles; ++i)
problem->particle_element_id[i] = -1;
int *old_particle_element_id = malloc(sizeof(int)*old_n_particles);
for (int i = 0; i < old_n_particles; i++){
old_particle_element_id[i] = -1;
}
double *old_particle_uvw = malloc(sizeof(double)*D*old_n_particles);
mesh_tree_particle_to_mesh(get_mesh_tree(problem), old_n_particles, old_particle_position, old_particle_element_id, old_particle_uvw);
compute_porosity(problem->mesh, problem->node_volume, problem->oldporosity, old_n_particles, old_particle_position, old_particle_volume, old_particle_element_id, old_particle_uvw, NULL);
mesh_tree_particle_to_mesh(get_mesh_tree(problem), problem->n_particles, problem->particle_position, problem->particle_element_id, problem->particle_uvw);
compute_porosity(problem->mesh, problem->node_volume, problem->porosity, problem->n_particles, problem->particle_position, problem->particle_volume, problem->particle_element_id, problem->particle_uvw, NULL);
free(old_particle_element_id);
free(old_particle_uvw);
}
#endif
void fluid_problem_set_coordinates(FluidProblem *problem, double *x) {
for (int i = 0; i < problem->mesh->n_nodes*3; ++i) {
problem->mesh->x[i] = x[i];
......
......@@ -320,7 +320,7 @@ class FluidProblem :
self.strong_cb_ref.append(bndcb)
self._lib.fluid_problem_set_strong_boundary(self._fp, tag.encode(), c_int(field), c_int(row), bndcb)
def adapt_mesh(self, lcmax, lcmin, n_el, old_n_particles, old_particle_position, old_particle_volume, cmax=1, cmin=1) :
def adapt_mesh(self, lcmax, lcmin, n_el, old_n_particles=0, old_particle_position=None, old_particle_volume=None, cmax=1, cmin=1) :
"""Adapts the mesh.
Keyword arguments:
......@@ -505,6 +505,7 @@ class FluidProblem :
self.solution()[:] -= sys.solve(rhs)
# print(time()-timer)
if check_residual_norm > 0:
# self._lib.fluid_problem_reset_boundary_force(self._fp)
self._lib.fluid_problem_assemble_system(self._fp,_np2c(solution_old),c_double(dt),_np2c(localv),_np2c(localm))
rhs = sys.local_to_global(localv,localm,self.solution().reshape([-1]),constraint_value)
norm = np.linalg.norm(rhs)
......
......@@ -7,42 +7,43 @@ from migflow import time_integration
from migflow import free_surface
outputdir = "output_2750"
outputdir = "output"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
outf = 2
outf = 40
h = 0.03
h = 0.01
l = 0.03
r = 0.001
xc = np.array([0.0,0])
# xc = np.array([0,-h/3+r])
vc = np.array([0.0,0.0])
xc = np.array([l/2, 1.5*r])
# xc = np.array([l/2, 1.75*r])
g = np.array([0, -9.81])
rho = 1200
rho_b = 2750#1200#2500#2750
rho_b = 2000#1200#2500#2750
m_b = rho_b*np.pi*r**2
print(f"m : {m_b} ---- r = {r}")
nu = 1.85e-5
dt = 1e-3
tEnd = 5
dt = 1e-4
tEnd = 1
fluid = fluid.FluidProblem(2,g,nu*rho,rho, order=[2,2,1], solver="petsc4py", solver_options="-pc_type lu")
fluid.load_msh("mesh.msh")
v_max = 0.1
v_max = 0.02 #0.01
def v_bubble(x, axis=0):
r = x[:,:]-xc[None,:]
r = r[:,:]/np.linalg.norm(r,axis=1)[:,None]
return v_max*r[:,axis]
# r = x[:,axis]-xc[None,axis]
# r /= np.linalg.norm(r,axis=0)
# return 0.1*r
return -v_max if axis==1 else 0
# r = x[:,:]-xc[None,:]
# r = r[:,:]/np.linalg.norm(r,axis=1)[:,None]
# return v_max*r[:,axis]
fluid.set_open_boundary("Lateral", velocity = [0,0])
fluid.set_open_boundary("Top", pressure=0)
......@@ -51,30 +52,49 @@ fluid.set_open_boundary("Bubble", velocity=[lambda x: v_bubble(x,0), lambda x :
fluid.set_strong_boundary("Bubble", 0, lambda x : v_bubble(x,0))
fluid.set_strong_boundary("Bubble", 1, lambda x : v_bubble(x,1))
nodes_bubble= np.unique(fluid.mesh_boundaries()[b"Bubble"])
nodes_boundary = np.unique(np.concatenate(list(fluid.mesh_boundaries().values()), axis=0))
x = fluid.coordinates()
t = 0
ii = 0
"""
def write_points():
with open("adapt/p.geo","w") as f:
f.write(f"Cx = {xc[0]};\n")
f.write(f"Cy = {xc[1]};\n")
f.write(f"Point(5) = {{{xc[0]},{xc[1]},0}};\n")
f.write(f"Point(6) = {{{xc[0]-r},{xc[1]},0,lc_circle}};\n")
f.write(f"Point(7) = {{{xc[0]+r},{xc[1]},0,lc_circle}};\n")
f.write(f"Point(8) = {{{xc[0]},{xc[1]-r},0,lc_circle}};\n")
f.write(f"Point(9) = {{{xc[0]},{xc[1]+r},0,lc_circle}};\n")
"""
nodes_bubble= np.unique(fluid.mesh_boundaries()[b"Bubble"])
nodes_boundary = np.unique(np.concatenate(list(fluid.mesh_boundaries().values()), axis=0))
fluid.write_vtk(outputdir, t, ii)
tic = time.time()
while(t < tEnd):
fluid.implicit_euler(dt, 1e-1)
if t>0.01:
xold = fluid.coordinates().copy()
fluid.implicit_euler(dt)
if ii > 10:
fbubble = fluid.boundary_force()[3,:] + m_b*g
print("\t fb = ", fluid.boundary_force()[3,1], "-----", m_b*g[1])
print("\t fbubble = ", fbubble)
vc = vc + fbubble/m_b*dt
vc[1] = vc[1] + fbubble[1]/m_b*dt
dxc = vc*dt
xold = x.copy()
R = min(0.75*(h/2-xc[1]), 0.75*(h/2+xc[1]))
free_surface.regularize_point(x, R, xc, r, dxc, nodes_boundary)
xc = xc+dxc
x = fluid.coordinates()
free_surface.regularize_point(x, xc[1], xc, r, dxc, nodes_boundary)
x[nodes_bubble,:2] = x[nodes_bubble,:2] + dxc[None,:2]
# xc[1] = xc[1] + dxc[1]
xc = xc + dxc
fluid.mesh_velocity()[:,:2] = (x[:,:2]-xold[:,:2])/dt
print("\t fb = ", fluid.boundary_force()[3,1], "-----", m_b*g[1], "-----------", fbubble)
print("\t v : ", vc)
t += dt
if ii%outf == 0:
ioutput = int(ii/outf)+1
......
H = 0.03;
H = 0.01;
L = 0.03;
r = 0.001;
Cx = 0.0;
// Cy = -H/3+r;
Cy = 0;
Cx = L/2;
Cy = 1.5*r;//1.5*r;
lc = 0.004;
lc_circle = 0.0005;
lc_circle= 0.0005;
Point(1) = {-L/2,-H/2,0,lc};
Point(2) = {L/2,-H/2,0,lc};
Point(3) = {L/2,H/2,0,lc};
Point(4) = {-L/2,H/2,0,lc};
Point(1) = {0,0,0,lc_circle};
Point(2) = {L,0,0,lc_circle};
Point(3) = {L,H,0,lc};
Point(4) = {0,H,0,lc};
Line(1) = {1,2};
Line(2) = {2,3};
Line(3) = {3,4};
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment