Skip to content
GitLab
Commit 8b4abd2e authored by dazerki's avatar dazerki
Browse files

link xmesh migflow

parent 2bb165b5
Pipeline #9821 passed with stages
in 3 minutes and 38 seconds
Hide whitespace changes
Inline Side-by-side
fluid/fluid_problem.c
View file @ 8b4abd2e
......@@ -246,7 +246,7 @@ static void f_boundary(WeakBoundary *wbnd, FluidProblem *problem,const double *n
if (unbnd - unmesh < 0) {
for (int id = 0; id < D; ++id) {
f0[U+id] += ((unbnd - unmesh)*(vid<0?0:uext[id])-(unold - unmesh)*u[id])*rho/c;
f00[(U+id)*n_fields+U+id] -= unold*rho/c;
f00[(U+id)*n_fields+U+id] -= (unold-unmesh)*rho/c;
}
}
}
......@@ -1336,6 +1336,7 @@ FluidProblem *fluid_problem_new(double *g, int n_fluids, double *mu, double *rho
problem->solution = NULL;
problem->mesh_velocity = NULL;
problem->concentration = NULL;
problem->concentration_cg = NULL;
problem->n_particles = 0;
problem->particle_uvw = NULL;
problem->particle_element_id = NULL;
......@@ -1366,6 +1367,7 @@ void fluid_problem_free(FluidProblem *problem)
free(problem->solution);
free(problem->mesh_velocity);
free(problem->concentration);
free(problem->concentration_cg);
free(problem->porosity);
free(problem->bulk_force);
free(problem->oldporosity);
......@@ -1662,6 +1664,8 @@ void fluid_problem_set_mesh(FluidProblem *problem, Mesh *mesh) {
if (problem->n_fluids == 2) {
free(problem->concentration);
problem->concentration = (double*)malloc(mesh->n_elements*N_N*sizeof(double));
free(problem->concentration_cg);
problem->concentration_cg = (double*)malloc(mesh->n_nodes*sizeof(double));
}
free(problem->mesh_velocity);
problem->mesh_velocity = (double*)calloc(mesh->n_nodes*D,sizeof(double));
......@@ -1904,6 +1908,12 @@ double *fluid_problem_concentration_dg_grad(const FluidProblem *p)
return p->grad_a_cg;
}
double *fluid_problem_concentration_cg(const FluidProblem *p)
{
fluid_problem_dg_to_cg(p, p->concentration, p->concentration_cg);
return p->concentration_cg;
}
double *fluid_problem_coordinates(const FluidProblem *p)
{
return p->mesh->x;
......
fluid/fluid_problem.h
View file @ 8b4abd2e
......@@ -68,6 +68,7 @@ struct FluidProblem {
double *porosity;
double *oldporosity;
double *concentration;
double *concentration_cg;
double *solution;
double *mesh_velocity;
double *node_volume;
......@@ -124,6 +125,7 @@ void fluid_problem_set_strong_boundary(FluidProblem *p, const char *tag, int fie
double *fluid_problem_old_porosity(const FluidProblem *p);
double *fluid_problem_porosity(const FluidProblem *p);
double *fluid_problem_concentration_dg(const FluidProblem *p);
double *fluid_problem_concentration_cg(const FluidProblem *p);
double *fluid_problem_element_size(const FluidProblem *p);
void fluid_problem_advance_concentration(FluidProblem *problem, double dt);
......
python/fluid.py
View file @ 8b4abd2e
......@@ -353,6 +353,7 @@ class FluidProblem :
("parent_node_id", self.parent_nodes())
]
if self._n_fluids == 2 :
data.append(("concentration_cg", self.concentration_cg()))
cell_data.append(("curvature",self.curvature()))
cell_data.append(("concentration",self.concentration_dg()))
cell_data.append(("stab",self._get_matrix("stab_param",self.n_elements(),1)))
......@@ -602,6 +603,10 @@ class FluidProblem :
"""Returns the concentration at discontinuous nodes"""
return self._get_matrix("concentration_dg", self.n_elements(), self._dim+1)
def concentration_cg(self):
"""Returns the concentration at discontinuous nodes"""
return self._get_matrix("concentration_cg", self.n_nodes(), 1)
def concentration_dg_grad(self):
"""Returns the concentration at discontinuous nodes"""
return self._get_matrix("concentration_dg_grad", self.n_nodes(), self._dim)
......
testcases/xmesh/bubble.py
View file @ 8b4abd2e
......@@ -8,6 +8,10 @@
# see <http://www.gnu.org/licenses/>.
#!/usr/bin/env python
# import sys
# sys.path.append('/home/Documents/migflow/migflow/python')
# import python
from migflow import fluid as mbfluid
from migflow import time_integration
......@@ -33,10 +37,9 @@ class Poiseuille(unittest.TestCase) :
t = 0
ii = 0
#physical parameters
g = np.array([0,-9.81]) # gravity
rho0 = 1 # fluid density
rho0 = 1 # fluid density
rho1 = 1000
nu0 = 1e-3 # kinematic viscosity
nu1 = 1e-3 # kinematic viscosity
......@@ -50,7 +53,7 @@ class Poiseuille(unittest.TestCase) :
#numerical parameters
lcmin = .1 # mesh size
dt = 0.005 # time step
dt = 0.0001 # time step
alpha = 1e-4 # stabilization coefficient
shutil.copy("mesh.msh", outputdir +"/mesh.msh")
......@@ -61,9 +64,10 @@ class Poiseuille(unittest.TestCase) :
fluid = mbfluid.FluidProblem(2,g,[nu0*rho0,nu1*rho1],[rho0,rho1])
fluid.load_msh("mesh.msh")
fluid.set_open_boundary("BottomMiddle", velocity=[0, boundaryV], concentration=1) #air entrance
fluid.set_open_boundary("BottomLeft",velocity=[0,0])
fluid.set_open_boundary("BottomRight",velocity=[0,0])
fluid.set_open_boundary("BottomMiddle", velocity=[0, boundaryV], concentration=1)
fluid.set_wall_boundary("BottomLeft",velocity=[0,0])
fluid.set_wall_boundary("BottomRight",velocity=[0,0])
# fluid.set_open_boundary("Insert", velocity=[0, boundaryV], concentration=1)
fluid.set_open_boundary("Top",pressure=0, concentration=0)
fluid.set_wall_boundary("Right",velocity=[0,0])
fluid.set_wall_boundary("Left",velocity=[0,0])
......@@ -79,28 +83,23 @@ class Poiseuille(unittest.TestCase) :
mesh,x = xmesh2d.Mesh.import_from_file("mesh.msh")
tracker = xmesh2d.FrontTracker(mesh,x)
ii = 0
tic = time.time()
while ii < 20:
while ii < 1000:
#Fluid solver
time_integration.iterate(fluid,None,dt)
t += dt
# level = fluid.concentration_dg_grad() - 0.5
# print(fluid.n_nodes())
# print(level)
concentration = fluid.concentration_dg_grad()[:,0]
# concentration = concentration.reshape((fluid.n_nodes(),1))
concentration = fluid.concentration_cg().reshape((fluid.n_nodes(),))
level = concentration - 0.5;
# print(level)
x_old = np.copy(x)
x = tracker.move_front(x, level)
x = tracker.relax(x, dt*5)
print(np.sum(x-x_old))
fluid.mesh_velocity()[:] = (x[:,:2] - fluid.coordinates()[:,:2])/dt
fluid.coordinates()[:] = x[:]
newx = np.copy(fluid.coordinates()[:])
newx = tracker.move_front(newx, level[:,])
newx = tracker.relax(newx, dt*5)
print(x[:,:2].shape)
print(np.sum(newx-fluid.coordinates()))
fluid.mesh_velocity()[:] = (newx[:,:2] - fluid.coordinates()[:,:2])/dt
fluid.coordinates()[:] = newx[:]
#Output files writting
if ii %outf == 0 :
......
testcases/xmesh/xmesh2d.py
View file @ 8b4abd2e
......@@ -6,7 +6,7 @@ gmsh.initialize()
class TrackerUI:
def __init__(self, initial_dt=0.05,pause=False):
def __init__(self, initial_dt=0.05):
gmsh.option.set_number("General.Antialiasing",1)
gmsh.option.set_number("General.SmallAxes",0)
gmsh.option.set_number("Print.Background",1)
......@@ -18,9 +18,8 @@ class TrackerUI:
gmsh.option.set_number(f"View[{self.view}].PointType",1)
gmsh.option.set_number(f"View[{self.view}].PointSize",5)
gmsh.option.set_number(f"View[{self.view}].ColormapBias",-0.15)
pause = 1 if pause else 0
gmsh.onelab.set(f"""[
{{ "type":"number", "name":"Pause", "values":[{pause}], "choices":[0, 1] }},
{{ "type":"number", "name":"Pause", "values":[0], "choices":[0, 1] }},
{{ "type":"number", "name":"Record", "values":[0], "choices":[0, 1] }},
{{ "type":"number", "name":"dt", "values":[{initial_dt}]}}
]""")
......@@ -28,11 +27,6 @@ class TrackerUI:
self.dt = gmsh.onelab.get_number("dt")[0]
gmsh.fltk.initialize()
def process_events(self):
gmsh.graphics.draw()
while gmsh.onelab.get_number("Pause") == 1:
gmsh.graphics.draw()
def update(self, x, front):
for tn in range(x.shape[0]):
gmsh.model.mesh.set_node(tn+1, x[tn], [0,0])
......@@ -166,6 +160,11 @@ class FrontTracker():
target = neighbours[np.argmin(sign[tn]*grad)]
alpha = np.clip((levelset[tn])/(levelset[tn]-levelset[target]),0,1)
xn[tn] = alpha*x[target]+(1-alpha)*x[tn]
# if alpha>0.95:
# print("HERE ")
# return True
# else:
# return False
return alpha == 1
else:
return False
......
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