Commit 45c2f7bc authored by Matthieu Constant's avatar Matthieu Constant
Browse files

change temporal for dt=0

parent 59cdf022
Pipeline #7758 passed with stages
in 5 minutes and 53 seconds
......@@ -284,8 +284,8 @@ void fluid_problem_compute_stab_parameters(FluidProblem *problem, double dt) {
fluid_problem_interpolate_rho_and_nu(problem,a, &rho, &mu);
double nu = mu/rho;
problem->taup[iel] = 1/sqrt(
(problem->temporal?pow2(2/dt):0.)
+(problem->advection ?pow2(2*normu/h):0.)
(dt!=0?pow2(2/dt):0.)
+((problem->advection & dt!=0 )?pow2(2*normu/h):0.)
+pow2(4*nu/pow2(h)));
problem->tauc[iel] = problem->advection?h*normu*fmin(h*normu/(6*nu),0.5):0.;
}
......@@ -390,7 +390,7 @@ static void fluid_problem_f(const FluidProblem *problem, const double *sol, doub
divu += du[iD][iD];
}
nold = sqrt(nold);
f0[P] = (c-cold)/dt;// + (sol[P]-solold[P])/dt*0.1;
if (dt!=0) {f0[P] = (c-cold)/dt;}// + (sol[P]-solold[P])/dt*0.1;
//f00[P*n_fields+P] = 1/dt*0.1;
double *g = problem->g;
double rhoreduced = (problem->reduced_gravity ? (rho-problem->rho[0]) : rho);
......@@ -403,12 +403,12 @@ static void fluid_problem_f(const FluidProblem *problem, const double *sol, doub
+ drag*u[i];
//+ 5.3e5*u[i];
f00[(U+i)*n_fields+U+i] = drag;//5.3e5;
if(problem->temporal) {
if(dt!=0) {
f0[U+i] += rho*(u[i]-uold[i])/dt;
f00[(U+i)*n_fields+U+i] += rho/dt;
}
f01[(U+i)*n_fields*D+P*D+i] = c;
if (problem->advection) {
if (problem->advection & dt!=0) {
// advection :
// dj(uj ui /c) = uj/c dj(ui) + ui/c dj(uj) - uj ui /(c*c) dj(c)
// = ui/c * (dj(uj)- uj/c*dj(c)) + uj/c*dj(ui)
......@@ -433,7 +433,7 @@ static void fluid_problem_f(const FluidProblem *problem, const double *sol, doub
f11[((U+i)*D+j)*n_fields*D+(U+i)*D+k] += supg*f01[(U+i)*n_fields*D+(U+i)*D+k];
}
double lsic = tauc*rho;
f1[(U+i)*D+i] += ((c-cold)/dt+divu)*lsic;
f1[(U+i)*D+i] += ((dt!=0?(c-cold)/dt:0)+divu)*lsic;
for (int j = 0; j < D; ++j) {
f11[((U+i)*D+i)*(n_fields*D)+(U+j)*D+j] += lsic;
}
......@@ -1248,7 +1248,7 @@ static void compute_porosity(Mesh *mesh, double *node_volume, double *porosity,
}
}
FluidProblem *fluid_problem_new(double *g, int n_fluids, double *mu, double *rho, double sigma, double volume_drag, double quadratic_drag, int drag_in_stab, double drag_coeff_factor, double ip_factor, int temporal, int advection)
FluidProblem *fluid_problem_new(double *g, int n_fluids, double *mu, double *rho, double sigma, double volume_drag, double quadratic_drag, int drag_in_stab, double drag_coeff_factor, double ip_factor, int advection)
{
if (n_fluids != 1 && n_fluids != 2) {
printf("error : only 1 or 2 fluids are supported\n");
......@@ -1263,7 +1263,6 @@ FluidProblem *fluid_problem_new(double *g, int n_fluids, double *mu, double *rho
problem->reduced_gravity = 0;
problem->strong_boundaries = NULL;
problem->n_strong_boundaries = 0;
problem->temporal = temporal;
problem->bulk_force = NULL;
problem->advection = advection;
problem->weak_boundaries = NULL;
......
......@@ -97,7 +97,6 @@ struct FluidProblem {
int reduced_gravity;
double stab_param;
int advection;
int temporal;
int drag_in_stab;
};
......@@ -107,7 +106,7 @@ void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, dou
void fluid_problem_set_bulk_force(FluidProblem *problem, double *force);
void fluid_problem_set_concentration_cg(FluidProblem *problem, double *concentration);
void fluid_problem_free(FluidProblem *problem);
FluidProblem *fluid_problem_new(double *g, int n_fluids, double *mu, double *rho, double sigma, double volume_drag, double quadratic_drag, int drag_in_stab, double drag_coeff_factor, double ip_factor,int temporal, int advection);
FluidProblem *fluid_problem_new(double *g, int n_fluids, double *mu, double *rho, double sigma, double volume_drag, double quadratic_drag, int drag_in_stab, double drag_coeff_factor, double ip_factor, int advection);
void fluid_problem_load_msh(FluidProblem *problem, const char *mesh_file_name);
void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin, double n_el, double cmax, double cmin, int old_n_particles, double *old_particle_position, double *old_particle_volume);
int *fluid_problem_particle_element_id(FluidProblem *problem);
......
......@@ -77,7 +77,7 @@ def _np2c(a,dtype=float) :
class FluidProblem :
"""Create numerical representation of the fluid."""
def __init__(self, dim, g, mu, rho, sigma=0, volume_drag=0., quadratic_drag=0., petsc_solver_type="", drag_in_stab=1, drag_coefficient_factor=1, ip_factor=1, temporal=True, advection=True):
def __init__(self, dim, g, mu, rho, sigma=0, volume_drag=0., quadratic_drag=0., petsc_solver_type="", drag_in_stab=1, drag_coefficient_factor=1, ip_factor=1, advection=True):
"""Build the fluid structure.
Keyword arguments:
......@@ -91,7 +91,6 @@ class FluidProblem :
coeff_stab -- optional argument used as coefficient for extra diffusivity added to stabilise the advection of concentration (only for two fluid problem)
drag_in_stab -- state if the drag force is in the stabilisation term
drag_coefficient_factor -- factor multiplicating the drag coefficient
temporal -- enable d/dt (i.e. False = steady)
advection -- enable advective terms (i.e. False = Stokes, True = Navier-Stokes)
......@@ -113,7 +112,7 @@ class FluidProblem :
n_fluids = np.require(rho,"float","C").reshape([-1]).shape[0]
self._n_fluids = n_fluids
self._dim = dim
self._fp = c_void_p(self._lib.fluid_problem_new(_np2c(g), n_fluids, _np2c(mu), _np2c(rho), c_double(sigma), c_double(volume_drag), c_double(quadratic_drag), c_int(drag_in_stab),c_double(drag_coefficient_factor),c_double(ip_factor),c_int(temporal), c_int(advection)))
self._fp = c_void_p(self._lib.fluid_problem_new(_np2c(g), n_fluids, _np2c(mu), _np2c(rho), c_double(sigma), c_double(volume_drag), c_double(quadratic_drag), c_int(drag_in_stab),c_double(drag_coefficient_factor),c_double(ip_factor), c_int(advection)))
if self._fp == None :
raise NameError("Cannot create fluid problem.")
......
......@@ -138,16 +138,23 @@ def iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e-8, external_particl
raise ValueError("fluid and particles are both None: not possible to iterate!")
if particles is not None and external_particles_forces is not None:
if external_particles_forces.shape!=particles.velocity().shape:
raise ValueError("external_particles_forces must have shape (number of particles,dimension!")
raise ValueError("external_particles_forces must have shape (number of particles,dimension)!")
if (particles is None or particles.r() is None) and fluid is not None:
fluid.implicit_euler(dt)
# For two fluids flows, advance the concentration field with the fluid velocity.
# The number of sub-time steps for the advection is automatically computed.
if fluid.n_fluids() == 2 :
if dt == 0:
raise ValueError("Time-step cannot be zero for two-fluids problems!")
fluid.advance_concentration(dt)
if particles is not None and fluid is None:
if dt == 0:
raise ValueError("Time-step cannot be zero for problems using grains!")
_advance_particles(particles,external_particles_forces,dt,min_nsub,contact_tol,after_sub_iter=after_sub_iter,max_split=max_split)
if fluid is not None and particles is not None and particles.r() is not None:
if dt == 0:
raise ValueError("Time-step cannot be zero for problems using grains!")
if fixed_grains:
fluid.implicit_euler(dt)
fluid.set_particles(particles.mass(), particles.volume(), particles.position(), particles.velocity(),-fluid.compute_node_force(dt))
......
L = 1;
L = .5;
H = 1;
y = 0;
lc = 0.1;
lc = .05;
Point(1) = {-L, H, 0,lc};
Point(2) = {-L, -H, 0,lc};
Point(3) = {L, -H, 0,lc};
Point(4) = {L, H, 0,lc};
Point(1) = {-L, H, 0, lc};
Point(2) = {-L, 0, 0, lc};
Point(3) = {L, 0, 0, lc};
Point(4) = {L, H, 0, lc};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
L1 = 1.5;
H1 = 3;
y = 0;
lc = 0.1;
Point(5) = {-L1, H1, 0,lc};
Point(6) = {-L1, -H1, 0,lc};
Point(7) = {L1, -H1, 0,lc};
Point(8) = {L1, H1, 0,lc};
Line(5) = {5, 6};
Line(6) = {6, 7};
Line(7) = {7, 8};
Line(8) = {8, 5};
Line Loop(1) = {1:4};
Line Loop(2) = {5:8};
Plane Surface(1) = {2,1};
Physical Line("BottomIn") = {2};
Physical Line("LeftIn") = {1};
Physical Line("RightIn") = {3};
Physical Line("TopIn") = {4};
Physical Line("BottomExt") = {6};
Physical Line("LeftExt") = {5};
Physical Line("RightExt") = {7};
Physical Line("TopExt") = {8};
Plane Surface(1) = {1};
Physical Line("Left") = {1};
Physical Line("Right") = {3};
Physical Line("Bottom") = {2};
Physical Line("Top") = {4};
Physical Surface("Domain") = {1};
Physical Point("PtFix") = {2};
Mesh.MshFileVersion = 2;
......@@ -21,7 +21,6 @@
#!/usr/bin/env python
from migflow import fluid as mbfluid
from migflow import scontact
from migflow import time_integration
import numpy as np
......@@ -32,72 +31,32 @@ import shutil
import random
def genInitialPosition(filename, r, L, L1, N, rhop) :
p = scontact.ParticleProblem(2)
#Loading of the mesh.msh file specifying physical boundaries name
p.load_msh_boundaries("mesh.msh", ["TopIn", "LeftIn", "RightIn","BottomIn","TopExt", "LeftExt", "RightExt","BottomExt"])
#Definition of the points where the grains are located
e = L/N**.5
#if e<2*r:
# print("Compacity too high!")
# exit(0)
x = np.arange(-L+r, L-r, e)
y = np.arange(-L1+r, -L-r, e)
x, y = np.meshgrid(x, y)
x = x.flat
y = y.flat
for i in range(len(x)) :
p.add_particle((x[i], y[i]), r, r**2 * np.pi * rhop);
print(p.position())
p.write_vtk(filename,0,0)
# 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)
j = 0
#physical parameters
g = np.array([0,0]) # gravity
rho = 1000 # fluid density
rhop = 1500 # grains density
L = 1 # half inner height
L1 = 3 # half outer height
mu = 1 # dynamic viscosity
r = 5e-3 # grains radius
#numerical parameters
compacity = .5 # compacity of the fixed granular pack
a_g = np.pi*r**2
a_b = L**2
N = max(compacity*a_b/a_g,1)
#Object particles creation
genInitialPosition(outputdir, r, L, L1, N, rhop)
p = scontact.ParticleProblem(2)
p.read_vtk(outputdir,0)
ii = 0
fluid = mbfluid.FluidProblem(2,g,mu,rho,temporal=False,advection=False)
vUp = 1
#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)
fluid = mbfluid.FluidProblem(2,g,mu,rho,advection=False)
fluid.load_msh("mesh.msh")
V = 0.001
fluid.set_wall_boundary("LeftIn",velocity=[0,-V])
fluid.set_wall_boundary("RightIn",velocity=[0,V])
fluid.set_wall_boundary("LeftExt",velocity=[0,0])
fluid.set_wall_boundary("RightExt",velocity=[0,0])
fluid.set_wall_boundary("BottomIn",velocity=[V,0])
fluid.set_wall_boundary("BottomExt",velocity=[0,0])
fluid.set_wall_boundary("TopExt",pressure=0)
fluid.set_wall_boundary("TopIn",velocity=[-V,0])
fp = np.zeros_like(fluid.porosity())
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces(),init=True)
fluid.export_vtk(outputdir,0,0)
t = 0
ii = 0
tic = time.time()
time_integration.iterate(fluid,p,1,fixed_grains=True)
t += 1
#Output files writting
p.write_vtk(outputdir, 1, 1)
fluid.export_vtk(outputdir, 1, 1)
ii += 1
fluid.set_wall_boundary("Left",velocity=[0,0])
fluid.set_wall_boundary("Bottom",velocity=[0,0])
fluid.set_wall_boundary("Top",pressure=0,velocity=[vUp,0])
fluid.set_wall_boundary("Right",velocity=[0,0])
fluid.export_vtk(outputdir,0,0)
#Fluid solver
time_integration.iterate(fluid,None,0)
#Output files writting
fluid.export_vtk(outputdir, 1, 1)
\ No newline at end of file
Markdown is supported
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