Commit 6093aafc authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

poiseuille fixed

parent 7bdeb284
......@@ -239,6 +239,42 @@ static void compute_node_force_implicit(FluidProblem *problem, double dt, double
typedef void f_cb(FluidProblem *problem, const double *n, double *f, const double *s, const double *ds, const double c, const double dt);
static void f_inflow(FluidProblem *problem,const double *n, double *f,const double *s,const double *ds,const double c,const double dt){
int n_fluids = problem->n_fluids;
int P = n_fluids*(D+1);
double p = s[P];
for (int ifluid = 0; ifluid < n_fluids; ++ifluid) {
int U = ifluid*(D+1);
int Q = ifluid*(D+1)+D;
double q = s[Q];
f[Q] = 0;
for (int id = 0; id < D; ++id) {
f[Q] -= s[U+id]*n[id];
//f[Q] += dp[id]*n[id]*epsilon*q;
f[U+id] = 0;//-p*q*n[id];
}
}
f[P] = 0;
}
static void f_outflow(FluidProblem *problem,const double *n, double *f,const double *s,const double *ds,const double c,const double dt){
int n_fluids = problem->n_fluids;
int P = n_fluids*(D+1);
double p = s[P];
for (int ifluid = 0; ifluid < n_fluids; ++ifluid) {
int U = ifluid*(D+1);
int Q = ifluid*(D+1)+D;
double q = s[Q];
f[Q] = 0;
for (int id = 0; id < D; ++id) {
f[Q] -= s[U+id]*n[id];
//f[Q] += dp[id]*n[id]*epsilon*q;
f[U+id] = 0;//-p*q*n[id];
}
}
f[P] = 0;
}
static void f_wall_pressure(FluidProblem *problem,const double *n, double *f,const double *s,const double *ds,const double c,const double dt)
{
int n_fluids = problem->n_fluids;
......@@ -255,6 +291,7 @@ static void f_wall_pressure(FluidProblem *problem,const double *n, double *f,con
double q = s[Q];
f[Q] = 0;
for (int id = 0; id < D; ++id) {
//f[Q] += dp[id]*n[id]*epsilon*q;
f[U+id] = q*p*n[id];
}
}
......@@ -301,7 +338,8 @@ static void fluid_problem_f(FluidProblem *problem, double *f0, double *f1, doubl
f0[Q] = (q-qold)/dt;
f0[P] -= q;
for (int i = 0; i < D; ++i) {
f0[U+i] = rho*((u[i]-uold[i])/dt + u[i]*oq2*divu + utau[i]*oq2) - p*dq[i] + (i==(D-1) ? -rho*problem->g : 0);
//f0[Q] += epsilonp*dq[i]*dp[i];
f0[U+i] = rho*((u[i]-uold[i])/dt + u[i]*oq2*divu + utau[i]*oq2) -p*dq[i]+ (i==(D-1) ? -q*rho*problem->g : 0);
for (int j = 0; j < D; ++j) {
f1[(U+j)*D+i] = mu*(tau[i][j]+tau[j][i]) + (i==j ? -q*p : 0);
}
......@@ -333,6 +371,14 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, d
f = f_wall_pressure;
else if (strcmp(mesh->phys_name[iphys],"Top") == 0)
f = f_wall_pressure;
else if (strcmp(mesh->phys_name[iphys],"LeftUp") == 0)
f = f_inflow;
else if (strcmp(mesh->phys_name[iphys],"LeftDown") == 0)
f = f_inflow;
else if (strcmp(mesh->phys_name[iphys],"RightUp") == 0)
f = f_outflow;
else if (strcmp(mesh->phys_name[iphys],"RightDown") == 0)
f = f_outflow;
else
continue;
......@@ -525,7 +571,6 @@ static void fluid_problem_volume(FluidProblem *problem, const double *solution_o
}
}
}
}
}
......@@ -596,8 +641,12 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
}
}
hxtLinearSystemAddToMatrix(lsys,iel,iel,local_matrix);
/*for (int i = 0; i < local_size; ++i)
printf("%8.2g ", local_vector[i]);
printf("\n");*/
hxtLinearSystemAddToRhs(lsys,rhs,iel,local_vector);
}
//exit(0);
free(forced_row);
free(all_local_matrix);
free(all_local_vector);
......
# Marblesbag - Copyright (C) <2010-2018>
# <Universite catholique de Louvain (UCL), Belgium
# Universite de Montpellier, France>
#
# List of the contributors to the development of Marblesbag: see AUTHORS file.
# Description and complete License: see LICENSE file.
#
# This program (Marblesbag) is free software:
# you can redistribute it and/or modify it under the terms of the GNU Lesser General
# Public License as published by the Free Software Foundation, either version
# 3 of the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with this program (see COPYING and COPYING.LESSER files). If not,
# see <http://www.gnu.org/licenses/>.
#!/usr/bin/env python
from marblesbag import fluid as fluid
from marblesbag import scontact2
import numpy as np
import os
import time
import shutil
import random
#Physical parameters for the drops are the ones presented by Metzger et al. (2007) "Falling clouds of particles in viscous fluids"
outputdir = "output"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
t = 0
ii = 0
#physical parameters
g = -9.81 # gravity
rho = 1000 # fluid density
nu = 1e-3 # kinematic viscosity
mu = nu*rho # dynamic viscosity
tEnd = 100000 # final time
#numerical parameters
lcmin = .1 # mesh size
dt = 5 # time step
alpha = 1e-7 # stabilization coefficient
epsilon = alpha*lcmin**2 /nu # stabilization parametre
print('epsilon',epsilon)
shutil.copy("mesh.msh", outputdir +"/mesh.msh")
outf = 1 # number of iterations between output files
#Object fluid creation + Boundary condition of the fluid (field 0 is horizontal velocity; field 1 is vertical velocity; field 2 is pressure)
#Format: strong_boundaries = [(Boundary tag, Fluid field, Value)
strong_boundaries = [
("Top",0,0,1),
("PtFix",2,3,0),
("Bottom",0,0,0),
("Bottom",1,1,0),
("Left",0,0,0),
("Left",1,1,0),
("Right",0,0,0),
("Right",1,1,0)
]
fluid = fluid.fluid_problem("mesh.msh",g,[nu*rho],[rho],epsilon,strong_boundaries,1)
ii = 0
t = 0
#set initial_condition
x = fluid.coordinates()
s = fluid.solution()
s[:,2] = 1.
#s[:,3] = -rho*g*(1-x[:,1])
fluid.export_vtk(outputdir,0,0)
ii = 0
tic = time.clock()
while ii < 100 :
#Fluid solver
fluid.implicit_euler(dt)
t += dt
#Output files writting
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
fluid.export_vtk(outputdir, t, ioutput)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.clock() - tic))
s = fluid.solution()
x = fluid.coordinates()
vel = (s[:,0]-1/(20*nu*rho)*x[:,1]*(1-x[:, 1]))**2
print('Error', (vel.sum())**.5)
if (vel.sum())**.5<1.5e-3:
exit(1)
else:
exit(0)
# Marblesbag - Copyright (C) <2010-2018>
# <Universite catholique de Louvain (UCL), Belgium
# Universite de Montpellier, France>
#
# List of the contributors to the development of Marblesbag: see AUTHORS file.
# Description and complete License: see LICENSE file.
#
# This program (Marblesbag) is free software:
# you can redistribute it and/or modify it under the terms of the GNU Lesser General
# Public License as published by the Free Software Foundation, either version
# 3 of the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with this program (see COPYING and COPYING.LESSER files). If not,
# see <http://www.gnu.org/licenses/>.
#!/usr/bin/env python
from marblesbag import fluid as fluid
from marblesbag import scontact2
import numpy as np
import os
import time
import shutil
import random
#Physical parameters for the drops are the ones presented by Metzger et al. (2007) "Falling clouds of particles in viscous fluids"
outputdir = "output"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
t = 0
ii = 0
#physical parameters
g = -9.81 # gravity
rho0 = 1050 # fluid density
rho1 = 1000
nu0 = 1e-3 # kinematic viscosity
nu1 = 1e-3 # kinematic viscosity
tEnd = 100000 # final time
#numerical parameters
lcmin = .1 # mesh size
dt = 1 # time step
alpha = 1e-6 # stabilization coefficien
epsilon = alpha*lcmin**2 /min(nu0,nu1) # stabilization parametre
print('epsilon',epsilon)
shutil.copy("mesh.msh", outputdir +"/mesh.msh")
outf = 1 # number of iterations between output files
#Object fluid creation + Boundary condition of the fluid (field 0 is horizontal velocity; field 1 is vertical velocity; field 2 is pressure)
#Format: strong_boundaries = [(Boundary tag, Fluid field, Value)
u0,v0,q0 = 0,1,2
u1,v1,q1 = 3,4,5
p = 6
strong_boundaries = [
#("Bottom",u0,u0,0.),
#("Bottom",u1,u1,0.),
#("Bottom",v0,v0,0.),
#("Bottom",v1,v1,0.),
#("Top",u0,u0,0.),
#("Top",u1,u1,0.),
#("Top",v0,v0,0.),
#("Top",v1,v1,0.),
#("Bottom",q1,q1,1.),
#("Bottom",q0,q0,0.),
#("Top",q1,q1,0.),
#("Top",q0,q0,1.),
("Top",p,p,0.),
#("Left",u0,u0,0),
#("Left",u1,u1,0),
#("Left",v0,v0,0),
#("Left",v1,v1,0),
#("Right",u0,u0,0),
#("Right",u1,u1,0),
#("Right",v0,v0,0),
#("Right",v1,v1,0),
]
fluid = fluid.fluid_problem("mesh.msh",g,[nu0*rho0,nu1*rho1],[rho0,rho1],epsilon,strong_boundaries,2)
ii = 0
t = 0
#set initial_condition
x = fluid.coordinates()
s = fluid.solution()
s[:,2] = 0.5#x[:,1]
s[:,5] = 1-s[:,2]
s[:,6] = -(rho0+rho1)/2*g*(1-x[:,1])
fluid.export_vtk(outputdir,0,0)
ii = 0
tic = time.clock()
while t < tEnd :
#Fluid solver
fluid.implicit_euler(dt)
t += dt
s = fluid.solution()
x = fluid.coordinates()
vel = (s[:,0]+s[:,3]-1/(20*nu0*rho0)*x[:,1]*(1-x[:, 1]))**2
print('Error', (vel.sum())**.5)
#Output files writting
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
fluid.export_vtk(outputdir, t, ioutput)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.clock() - tic))
......@@ -59,9 +59,18 @@ outf = 1 # number of iterations between ou
#Object fluid creation + Boundary condition of the fluid (field 0 is horizontal velocity; field 1 is vertical velocity; field 2 is pressure)
#Format: strong_boundaries = [(Boundary tag, Fluid field, Value)
strong_boundaries = [("Top",1,1,0.),("Bottom",1,1,0.),
("Left",0,0,0),("Right",0,0,0),
("Top",2,3,0)]
strong_boundaries = [
#("Top",0,0,0.),
#("Top",1,1,0.),
#("Bottom",0,0,0.),
#("Bottom",1,1,0.),
#("Left",0,0,0),
#("Left",1,1,0),
#("Right",0,0,0),
#("Right",1,1,0),
("Top",2,3,0),
#("Top",2,2,1),
]
fluid = fluid.fluid_problem("mesh.msh",g,[nu*rho],[rho],epsilon,strong_boundaries,1)
ii = 0
t = 0
......@@ -70,7 +79,7 @@ t = 0
x = fluid.coordinates()
s = fluid.solution()
s[:,2] = 1.
#s[:,3] = -rho*g*(1-x[:,1])
s[:,3] = -rho*g*(1-x[:,1])
fluid.export_vtk(outputdir,0,0)
......
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