Skip to content
GitLab
Projects
Groups
Snippets
/
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
fluidparticles
MigFlow
Commits
8591a6ba
Commit
8591a6ba
authored
Nov 10, 2021
by
Jonathan Lambrechts
Browse files
usolid
parent
828e88d5
Pipeline
#9953
failed with stages
in 45 minutes and 39 seconds
Changes
6
Pipelines
2
Hide whitespace changes
Inline
Side-by-side
fluid/fluid_problem.c
View file @
8591a6ba
...
...
@@ -200,7 +200,7 @@ static void f_boundary(WeakBoundary *wbnd, FluidProblem *problem,const double *n
dp
[
iD
]
=
ds
[
P
*
D
+
iD
];
unold
+=
sold
[
U
+
iD
]
*
n
[
iD
];
unmesh
+=
mesh_velocity
[
iD
]
*
n
[
iD
];
uext
[
iD
]
=
(
vid
<
0
?
s
[
U
+
iD
]
:
data
[
vid
+
iD
]);
uext
[
iD
]
=
(
vid
<
0
?
s
[
U
+
iD
]
:
data
[
vid
+
iD
]
*
c
);
unext
+=
uext
[
iD
]
*
n
[
iD
];
dcn
+=
dc
[
iD
]
*
n
[
iD
];
s_c
+=
(
u
[
iD
]
-
uext
[
iD
])
*
n
[
iD
];
...
...
@@ -381,7 +381,7 @@ static void particle_force_f(FluidProblem *problem, double *f, double *dfdu, dou
*
dfddp
=
gammastar
*
dt
/
mass
*
vol
;
//-vol;
}
static
void
fluid_problem_f
(
const
FluidProblem
*
problem
,
const
double
*
sol
,
double
*
dsol
,
const
double
*
solold
,
const
double
*
dsolold
,
const
double
*
mesh_velocity
,
const
double
c
,
const
double
*
dc
,
const
double
*
bf
,
const
double
cold
,
const
double
rho
,
const
double
mu
,
double
dt
,
int
iel
,
double
*
f0
,
double
*
f1
,
double
*
f00
,
double
*
f10
,
double
*
f01
,
double
*
f11
)
{
static
void
fluid_problem_f
(
const
FluidProblem
*
problem
,
const
double
*
sol
,
double
*
dsol
,
const
double
*
solold
,
const
double
*
dsolold
,
const
double
*
mesh_velocity
,
const
double
c
,
const
double
*
dc
,
const
double
*
bf
,
const
double
cold
,
const
double
dus
[
D
][
D
],
const
double
rho
,
const
double
mu
,
double
dt
,
int
iel
,
double
*
f0
,
double
*
f1
,
double
*
f00
,
double
*
f10
,
double
*
f01
,
double
*
f11
)
{
double
p
=
sol
[
P
];
double
taup
=
problem
->
taup
[
iel
];
double
tauc
=
problem
->
tauc
[
iel
];
...
...
@@ -390,6 +390,7 @@ static void fluid_problem_f(const FluidProblem *problem, const double *sol, doub
if
(
problem
->
stab_param
!=
0
)
taup
=
tauc
=
0
;
double
divu
=
0
;
double
divus
=
0
;
double
nold
=
0
;
for
(
int
iD
=
0
;
iD
<
D
;
++
iD
)
{
u
[
iD
]
=
sol
[
U
+
iD
];
...
...
@@ -403,8 +404,15 @@ static void fluid_problem_f(const FluidProblem *problem, const double *sol, doub
}
divu
+=
du
[
iD
][
iD
];
}
if
(
problem
->
usolid
)
{
for
(
int
iD
=
0
;
iD
<
D
;
++
iD
)
divus
+=
dus
[
iD
][
iD
];
}
else
{
divus
=
(
c
-
cold
)
/
dt
;
}
nold
=
sqrt
(
nold
);
f0
[
P
]
=
(
c
-
cold
)
/
dt
;
f0
[
P
]
=
divus
;
double
*
g
=
problem
->
g
;
double
rhoreduced
=
(
problem
->
reduced_gravity
?
(
rho
-
problem
->
rho
[
0
])
:
rho
);
double
drag
=
problem
->
volume_drag
+
problem
->
quadratic_drag
*
nold
;
...
...
@@ -446,7 +454,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
]
+=
(
divus
+
divu
)
*
lsic
;
for
(
int
j
=
0
;
j
<
D
;
++
j
)
{
f11
[((
U
+
i
)
*
D
+
i
)
*
(
n_fields
*
D
)
+
(
U
+
j
)
*
D
+
j
]
+=
lsic
;
}
...
...
@@ -757,7 +765,7 @@ double fluid_problem_a_integ_bnd(FluidProblem *problem, double dt)
return
i_bnd
;
}
static
void
fluid_problem_volume
(
FluidProblem
*
problem
,
const
double
*
solution_old
,
double
dt
,
double
*
all_local_vector
,
double
*
all_local_matrix
)
static
void
fluid_problem_volume
(
FluidProblem
*
problem
,
const
double
*
solution_old
,
const
double
*
u_solid
,
double
dt
,
double
*
all_local_vector
,
double
*
all_local_matrix
)
{
const
Mesh
*
mesh
=
problem
->
mesh
;
const
double
*
porosity
=
problem
->
porosity
;
...
...
@@ -780,6 +788,7 @@ static void fluid_problem_volume(FluidProblem *problem, const double *solution_o
for
(
int
iel
=
0
;
iel
<
mesh
->
n_elements
;
++
iel
)
{
const
int
*
el
=
&
mesh
->
elements
[
iel
*
N_N
];
double
dxidx
[
D
][
D
],
dphi
[
N_N
][
D
];
double
dus
[
D
][
D
]
=
{
0
};
const
double
detj
=
mesh_dxidx
(
mesh
,
iel
,
dxidx
);
grad_shape_functions
(
dxidx
,
dphi
);
double
*
local_vector
=
&
all_local_vector
[
local_size
*
iel
];
...
...
@@ -802,6 +811,9 @@ static void fluid_problem_volume(FluidProblem *problem, const double *solution_o
}
for
(
int
k
=
0
;
k
<
D
;
++
k
){
dc
[
k
]
+=
problem
->
porosity
[
el
[
i
]]
*
dphi
[
i
][
k
];
for
(
int
j
=
0
;
j
<
D
;
++
j
)
{
dus
[
j
][
k
]
+=
u_solid
[
el
[
i
]
*
DIMENSION
+
j
]
*
dphi
[
i
][
k
];
}
}
if
(
problem
->
n_fluids
==
2
)
{
for
(
int
k
=
0
;
k
<
D
;
++
k
){
...
...
@@ -856,7 +868,7 @@ static void fluid_problem_volume(FluidProblem *problem, const double *solution_o
const
double
jw
=
QW
[
iqp
]
*
detj
;
problem
->
kinetic_energy
+=
rho
*
sqrt
(
s
[
U
]
*
s
[
U
]
+
s
[
U
+
1
]
*
s
[
U
+
1
])
*
sqrt
(
s
[
U
]
*
s
[
U
]
+
s
[
U
+
1
]
*
s
[
U
+
1
])
/
2
*
jw
;
fluid_problem_f
(
problem
,
s
,
ds
,
sold
,
dsold
,
mesh_velocity
,
c
,
dc
,
bf
,
cold
,
rho
,
mu
,
dt
,
iel
,
f0
,
f1
,
f00
,
f10
,
f01
,
f11
);
fluid_problem_f
(
problem
,
s
,
ds
,
sold
,
dsold
,
mesh_velocity
,
c
,
dc
,
bf
,
cold
,
dus
,
rho
,
mu
,
dt
,
iel
,
f0
,
f1
,
f00
,
f10
,
f01
,
f11
);
for
(
int
ifield
=
0
;
ifield
<
n_fields
;
++
ifield
)
{
for
(
int
iphi
=
0
;
iphi
<
N_SF
;
++
iphi
){
local_vector
[
iphi
+
N_SF
*
ifield
]
+=
phi
[
iphi
]
*
f0
[
ifield
]
*
jw
;
...
...
@@ -1169,77 +1181,46 @@ static void compute_porosity(Mesh *mesh, double *node_volume, double *porosity,
}
}
static
void
fluid_problem_update_porosity
(
FluidProblem
*
problem
,
double
dt
)
static
void
compute_u_solid
(
Mesh
*
mesh
,
double
*
node_volume
,
double
*
u_solid
,
int
n_particles
,
double
*
particle_volume
,
double
*
particle_velocity
,
int
*
particle_element_id
,
double
*
particle_uvw
)
{
double
*
old_p
=
malloc
(
sizeof
(
double
)
*
DIMENSION
*
problem
->
n_particles
);
for
(
int
i
=
0
;
i
<
problem
->
n_particles
*
DIMENSION
;
++
i
)
{
old_p
[
i
]
=
problem
->
particle_position
[
i
]
-
dt
*
problem
->
particle_velocity
[
i
];
}
double
*
old_uvw
=
malloc
(
sizeof
(
double
)
*
DIMENSION
*
problem
->
n_particles
);
mesh_tree_particle_to_mesh
(
get_mesh_tree
(
problem
),
problem
->
n_particles
,
old_p
,
problem
->
particle_element_id
,
old_uvw
);
compute_porosity
(
problem
->
mesh
,
problem
->
node_volume
,
problem
->
oldporosity
,
problem
->
n_particles
,
problem
->
particle_volume
,
problem
->
particle_element_id
,
old_uvw
);
int
*
old_eid
=
malloc
(
sizeof
(
int
)
*
problem
->
n_particles
);
for
(
int
i
=
0
;
i
<
problem
->
n_particles
;
++
i
)
{
old_eid
[
i
]
=
problem
->
particle_element_id
[
i
];
int
independent_nodes
=
mesh
->
n_nodes
-
mesh
->
n_periodic
;
double
*
periodic_u_solid
=
malloc
(
sizeof
(
double
)
*
independent_nodes
*
DIMENSION
);
for
(
int
i
=
0
;
i
<
independent_nodes
*
DIMENSION
;
++
i
){
periodic_u_solid
[
i
]
=
0
;
}
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_volume
,
problem
->
particle_element_id
,
problem
->
particle_uvw
);
double
*
mod_volume
=
malloc
(
sizeof
(
double
)
*
problem
->
n_particles
);
double
*
mod_porosity
=
malloc
(
sizeof
(
double
)
*
problem
->
mesh
->
n_nodes
);
int
*
mod_eid
=
malloc
(
sizeof
(
int
)
*
problem
->
n_particles
);
double
*
mod_uvw
=
malloc
(
sizeof
(
double
)
*
problem
->
n_particles
*
DIMENSION
);
int
mod_n
=
0
;
for
(
int
i
=
0
;
i
<
problem
->
n_particles
;
++
i
)
{
if
(
old_eid
[
i
]
<
0
&&
problem
->
particle_element_id
[
i
]
>=
0
)
{
mod_eid
[
mod_n
]
=
problem
->
particle_element_id
[
i
];
mod_volume
[
mod_n
]
=
problem
->
particle_volume
[
i
];
for
(
int
j
=
0
;
j
<
DIMENSION
;
++
j
)
{
mod_uvw
[
mod_n
*
DIMENSION
+
j
]
=
problem
->
particle_uvw
[
i
*
DIMENSION
+
j
];
}
mod_n
+=
1
;
}
else
if
(
old_eid
[
i
]
>=
0
&&
problem
->
particle_element_id
[
i
]
<
0
)
{
mod_eid
[
mod_n
]
=
old_eid
[
i
];
mod_volume
[
mod_n
]
=
-
problem
->
particle_volume
[
i
];
for
(
int
j
=
0
;
j
<
DIMENSION
;
++
j
)
{
mod_uvw
[
mod_n
*
DIMENSION
+
j
]
=
old_uvw
[
i
*
DIMENSION
+
j
];
int
*
periodic_mapping
=
mesh
->
periodic_mapping
;
for
(
int
i
=
0
;
i
<
n_particles
;
++
i
)
{
if
(
particle_element_id
[
i
]
==
-
1
)
continue
;
double
sf
[
N_SF
];
shape_functions
(
&
particle_uvw
[
i
*
D
],
sf
);
const
int
*
el
=
&
mesh
->
elements
[
particle_element_id
[
i
]
*
N_N
];
for
(
int
j
=
0
;
j
<
N_N
;
++
j
)
{
for
(
int
k
=
0
;
k
<
DIMENSION
;
++
k
)
{
periodic_u_solid
[
periodic_mapping
[
el
[
j
]]
*
DIMENSION
+
k
]
+=
sf
[
j
]
*
particle_volume
[
i
]
*
particle_velocity
[
i
*
DIMENSION
+
k
];
}
mod_n
+=
1
;
}
}
compute_porosity
(
problem
->
mesh
,
problem
->
node_volume
,
mod_porosity
,
mod_n
,
mod_volume
,
mod_eid
,
mod_uvw
);
int
n_fields
=
fluid_problem_n_fields
(
problem
);
for
(
int
i
=
0
;
i
<
problem
->
mesh
->
n_nodes
;
++
i
)
{
double
op
=
problem
->
oldporosity
[
i
];
double
np
=
op
-
1
+
mod_porosity
[
i
];
problem
->
oldporosity
[
i
]
=
np
;
for
(
int
d
=
0
;
d
<
DIMENSION
;
++
d
)
{
problem
->
solution
[
n_fields
*
i
+
U
+
d
]
*=
np
/
op
;
for
(
int
i
=
0
;
i
<
mesh
->
n_nodes
;
++
i
){
for
(
int
k
=
0
;
k
<
DIMENSION
;
++
k
)
{
if
(
node_volume
[
i
]
==
0
.){
u_solid
[
i
*
DIMENSION
+
k
]
=
0
.;
}
else
{
u_solid
[
i
*
DIMENSION
+
k
]
=
periodic_u_solid
[
periodic_mapping
[
i
]
*
DIMENSION
+
k
]
/
node_volume
[
i
];
}
}
}
free
(
old_eid
);
free
(
old_uvw
);
free
(
old_p
);
free
(
mod_volume
);
free
(
mod_porosity
);
free
(
mod_eid
);
free
(
mod_uvw
);
free
(
periodic_u_solid
);
}
void
fluid_problem_u_solid
(
FluidProblem
*
problem
,
double
*
u_solid
)
{
compute_u_solid
(
problem
->
mesh
,
problem
->
node_volume
,
u_solid
,
problem
->
n_particles
,
problem
->
particle_volume
,
problem
->
particle_velocity
,
problem
->
particle_element_id
,
problem
->
particle_uvw
);
}
void
fluid_problem_assemble_system
(
FluidProblem
*
problem
,
const
double
*
solution_old
,
double
dt
,
double
*
all_local_vector
,
double
*
all_local_matrix
)
{
const
Mesh
*
mesh
=
problem
->
mesh
;
fluid_problem_update_porosity
(
problem
,
dt
);
int
n_fields
=
fluid_problem_n_fields
(
problem
);
...
...
@@ -1270,7 +1251,10 @@ void fluid_problem_assemble_system(FluidProblem *problem, const double *solution
}
fluid_problem_node_force_volume
(
problem
,
solution_old
,
dt
,
all_local_vector
,
all_local_matrix
);
compute_weak_boundary_conditions
(
problem
,
solution_old
,
dt
,
all_local_vector
,
all_local_matrix
);
fluid_problem_volume
(
problem
,
solution_old
,
dt
,
all_local_vector
,
all_local_matrix
);
double
*
u_solid
=
malloc
(
sizeof
(
double
)
*
problem
->
mesh
->
n_nodes
*
DIMENSION
);
compute_u_solid
(
problem
->
mesh
,
problem
->
node_volume
,
u_solid
,
problem
->
n_particles
,
problem
->
particle_volume
,
problem
->
particle_velocity
,
problem
->
particle_element_id
,
problem
->
particle_uvw
);
fluid_problem_volume
(
problem
,
solution_old
,
u_solid
,
dt
,
all_local_vector
,
all_local_matrix
);
free
(
u_solid
);
if
(
problem
->
n_fluids
==
2
){
double
*
grad_a_cg
=
malloc
(
sizeof
(
double
)
*
mesh
->
n_nodes
*
D
);
fluid_problem_dg_to_cg_grad
(
problem
,
problem
->
concentration
,
problem
->
grad_a_cg
);
...
...
@@ -1356,7 +1340,7 @@ static void fluid_problem_compute_node_volume(FluidProblem *problem) {
free
(
periodic_node_volume
);
}
FluidProblem
*
fluid_problem_new
(
double
*
g
,
int
n_fluids
,
double
*
mu
,
double
*
rho
,
double
sigma
,
double
coeffStab
,
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
coeffStab
,
double
volume_drag
,
double
quadratic_drag
,
int
drag_in_stab
,
double
drag_coeff_factor
,
double
ip_factor
,
int
temporal
,
int
advection
,
int
usolid
)
{
if
(
n_fluids
!=
1
&&
n_fluids
!=
2
)
{
printf
(
"error : only 1 or 2 fluids are supported
\n
"
);
...
...
@@ -1368,6 +1352,7 @@ FluidProblem *fluid_problem_new(double *g, int n_fluids, double *mu, double *rho
problem
->
ip_factor
=
ip_factor
;
problem
->
stab_param
=
0
;
problem
->
drag_in_stab
=
drag_in_stab
;
problem
->
usolid
=
usolid
;
problem
->
reduced_gravity
=
0
;
problem
->
strong_boundaries
=
NULL
;
problem
->
n_strong_boundaries
=
0
;
...
...
@@ -1719,8 +1704,6 @@ void fluid_problem_set_mesh(FluidProblem *problem, Mesh *mesh) {
free
(
problem
->
boundary_force
);
problem
->
boundary_force
=
(
double
*
)
malloc
(
DIMENSION
*
mesh
->
n_boundaries
*
sizeof
(
double
));
for
(
int
i
=
0
;
i
<
mesh
->
n_boundaries
*
DIMENSION
;
i
++
)
problem
->
boundary_force
[
i
]
=
0
;
free
(
problem
->
solution
);
int
n_fields
=
fluid_problem_n_fields
(
problem
);
problem
->
solution
=
(
double
*
)
malloc
(
mesh
->
n_nodes
*
n_fields
*
sizeof
(
double
));
...
...
@@ -1813,10 +1796,6 @@ void fluid_problem_set_coordinates(FluidProblem *problem, double *x) {
}
void
fluid_problem_set_particles
(
FluidProblem
*
problem
,
int
n
,
double
*
mass
,
double
*
volume
,
double
*
position
,
double
*
velocity
,
double
*
contact
)
{
double
*
cold
=
malloc
(
sizeof
(
double
)
*
problem
->
mesh
->
n_nodes
);
compute_porosity
(
problem
->
mesh
,
problem
->
node_volume
,
cold
,
problem
->
n_particles
,
problem
->
particle_volume
,
problem
->
particle_element_id
,
problem
->
particle_uvw
);
if
(
problem
->
n_particles
!=
n
)
{
problem
->
n_particles
=
n
;
free
(
problem
->
particle_uvw
);
...
...
@@ -1846,19 +1825,20 @@ void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, dou
}
problem
->
particle_element_id
[
i
]
=-
1
;
}
double
*
cnew
=
malloc
(
sizeof
(
double
)
*
problem
->
mesh
->
n_nodes
);
mesh_tree_particle_to_mesh
(
get_mesh_tree
(
problem
),
n
,
position
,
problem
->
particle_element_id
,
problem
->
particle_uvw
);
compute_porosity
(
problem
->
mesh
,
problem
->
node_volume
,
cnew
,
double
*
c
=
malloc
(
sizeof
(
double
)
*
problem
->
mesh
->
n_nodes
);
mesh_tree_particle_to_mesh
(
get_mesh_tree
(
problem
),
problem
->
n_particles
,
position
,
problem
->
particle_element_id
,
problem
->
particle_uvw
);
compute_porosity
(
problem
->
mesh
,
problem
->
node_volume
,
c
,
problem
->
n_particles
,
problem
->
particle_volume
,
problem
->
particle_element_id
,
problem
->
particle_uvw
);
int
n_fields
=
fluid_problem_n_fields
(
problem
);
for
(
int
i
=
0
;
i
<
problem
->
mesh
->
n_nodes
;
++
i
)
{
for
(
int
d
=
0
;
d
<
DIMENSION
;
++
d
)
{
problem
->
solution
[
n_fields
*
i
+
U
+
d
]
*=
c
new
[
i
]
/
cold
[
i
];
problem
->
solution
[
n_fields
*
i
+
U
+
d
]
*=
c
[
i
]
/
problem
->
porosity
[
i
];
}
problem
->
porosity
[
i
]
=
c
[
i
];
}
free
(
cnew
);
free
(
cold
);
free
(
c
);
}
void
fluid_problem_move_particles
(
FluidProblem
*
problem
,
int
n
,
double
*
position
,
double
*
velocity
,
double
*
contact
)
...
...
@@ -1867,6 +1847,21 @@ void fluid_problem_move_particles(FluidProblem *problem, int n, double *position
printf
(
"the number of particle has changed and the
\"
set_particles
\"
function has not been called.
\n
"
);
exit
(
1
);
}
int
n_fields
=
fluid_problem_n_fields
(
problem
);
int
*
old_eid
=
malloc
(
sizeof
(
int
)
*
problem
->
n_particles
);
double
*
old_uvw
=
malloc
(
sizeof
(
double
)
*
problem
->
n_particles
*
DIMENSION
);
for
(
int
i
=
0
;
i
<
problem
->
n_particles
;
++
i
)
{
old_eid
[
i
]
=
problem
->
particle_element_id
[
i
];
for
(
int
j
=
0
;
j
<
DIMENSION
;
j
++
)
{
old_uvw
[
i
*
DIMENSION
+
j
]
=
problem
->
particle_uvw
[
i
*
DIMENSION
+
j
];
}
}
for
(
int
i
=
0
;
i
<
problem
->
mesh
->
n_nodes
;
++
i
)
{
problem
->
oldporosity
[
i
]
=
problem
->
porosity
[
i
];
}
for
(
int
i
=
0
;
i
<
n
;
++
i
)
{
for
(
int
k
=
0
;
k
<
D
;
++
k
)
{
problem
->
particle_position
[
i
*
D
+
k
]
=
position
[
i
*
D
+
k
];
...
...
@@ -1874,6 +1869,52 @@ void fluid_problem_move_particles(FluidProblem *problem, int n, double *position
problem
->
particle_contact
[
i
*
D
+
k
]
=
contact
[
i
*
D
+
k
];
}
}
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_volume
,
problem
->
particle_element_id
,
problem
->
particle_uvw
);
double
*
mod_volume
=
malloc
(
sizeof
(
double
)
*
problem
->
n_particles
);
double
*
mod_porosity
=
malloc
(
sizeof
(
double
)
*
problem
->
mesh
->
n_nodes
);
int
*
mod_eid
=
malloc
(
sizeof
(
int
)
*
problem
->
n_particles
);
double
*
mod_uvw
=
malloc
(
sizeof
(
double
)
*
problem
->
n_particles
*
DIMENSION
);
int
mod_n
=
0
;
for
(
int
i
=
0
;
i
<
problem
->
n_particles
;
++
i
)
{
if
(
old_eid
[
i
]
<
0
&&
problem
->
particle_element_id
[
i
]
>=
0
)
{
mod_eid
[
mod_n
]
=
problem
->
particle_element_id
[
i
];
mod_volume
[
mod_n
]
=
problem
->
particle_volume
[
i
];
for
(
int
j
=
0
;
j
<
DIMENSION
;
++
j
)
{
mod_uvw
[
mod_n
*
DIMENSION
+
j
]
=
problem
->
particle_uvw
[
i
*
DIMENSION
+
j
];
}
mod_n
+=
1
;
}
else
if
(
old_eid
[
i
]
>=
0
&&
problem
->
particle_element_id
[
i
]
<
0
)
{
mod_eid
[
mod_n
]
=
old_eid
[
i
];
mod_volume
[
mod_n
]
=
-
problem
->
particle_volume
[
i
];
for
(
int
j
=
0
;
j
<
DIMENSION
;
++
j
)
{
mod_uvw
[
mod_n
*
DIMENSION
+
j
]
=
old_uvw
[
i
*
DIMENSION
+
j
];
}
mod_n
+=
1
;
}
}
compute_porosity
(
problem
->
mesh
,
problem
->
node_volume
,
mod_porosity
,
mod_n
,
mod_volume
,
mod_eid
,
mod_uvw
);
for
(
int
i
=
0
;
i
<
problem
->
mesh
->
n_nodes
;
++
i
)
{
double
mod_compacity
=
(
1
-
mod_porosity
[
i
]);
double
f
=
(
problem
->
oldporosity
[
i
]
-
mod_compacity
)
/
problem
->
oldporosity
[
i
];
for
(
int
d
=
0
;
d
<
DIMENSION
;
++
d
)
{
problem
->
solution
[
n_fields
*
i
+
U
+
d
]
*=
f
;
}
problem
->
oldporosity
[
i
]
*=
f
;
}
free
(
old_eid
);
free
(
old_uvw
);
free
(
mod_volume
);
free
(
mod_porosity
);
free
(
mod_eid
);
free
(
mod_uvw
);
}
double
*
fluid_problem_particle_uvw
(
FluidProblem
*
problem
)
{
...
...
fluid/fluid_problem.h
View file @
8591a6ba
...
...
@@ -99,6 +99,7 @@ struct FluidProblem {
int
reduced_gravity
;
double
stab_param
;
int
advection
;
int
usolid
;
int
temporal
;
int
drag_in_stab
;
};
...
...
@@ -110,7 +111,7 @@ void fluid_problem_move_particles(FluidProblem *problem, int n, double *position
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
coeffStab
,
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
coeffStab
,
double
volume_drag
,
double
quadratic_drag
,
int
drag_in_stab
,
double
drag_coeff_factor
,
double
ip_factor
,
int
temporal
,
int
advection
,
int
usolid
);
void
fluid_problem_adapt_mesh
(
FluidProblem
*
problem
,
Mesh
*
new_mesh
,
int
old_n_particles
,
double
*
old_particle_position
,
double
*
old_particle_volume
);
void
fluid_problem_set_mesh
(
FluidProblem
*
problem
,
Mesh
*
mesh
);
int
*
fluid_problem_particle_element_id
(
FluidProblem
*
problem
);
...
...
python/fluid.py
View file @
8591a6ba
...
...
@@ -123,7 +123,7 @@ def _load_msh(mesh_file_name, lib, dim) :
class
FluidProblem
:
"""Creates the numerical representation of the fluid."""
def
__init__
(
self
,
dim
,
g
,
mu
,
rho
,
sigma
=
0
,
coeff_stab
=
0.01
,
volume_drag
=
0.
,
quadratic_drag
=
0.
,
solver
=
None
,
solver_options
=
""
,
drag_in_stab
=
1
,
drag_coefficient_factor
=
1
,
ip_factor
=
1
,
temporal
=
True
,
advection
=
True
):
def
__init__
(
self
,
dim
,
g
,
mu
,
rho
,
sigma
=
0
,
coeff_stab
=
0.01
,
volume_drag
=
0.
,
quadratic_drag
=
0.
,
solver
=
None
,
solver_options
=
""
,
drag_in_stab
=
1
,
drag_coefficient_factor
=
1
,
ip_factor
=
1
,
temporal
=
True
,
advection
=
True
,
usolid
=
False
):
"""Builds the fluid structure.
Keyword arguments:
...
...
@@ -138,7 +138,9 @@ class FluidProblem :
drag_in_stab -- States if the drag force is in the stabilisation term
drag_coefficient_factor -- Factor multiplicating the drag coefficient
temporal -- Enables d/dt (i.e. False = steady)
advection -- Enables advective terms (i.e. False = Stokes, True = Navier-Stokes)
advection -- Enables advective terms (i.e. False = Stokes, True = Navier-Stokes
usolid -- if False, div.u_solid is evaluated = (porosity-old_porosity)/dt
)
Raises:
...
...
@@ -161,7 +163,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
(
coeff_stab
),
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
(
coeff_stab
),
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
)
,
c_int
(
usolid
)
))
if
self
.
_fp
==
None
:
raise
NameError
(
"Cannot create fluid problem."
)
...
...
@@ -341,11 +343,14 @@ class FluidProblem :
if
self
.
_dim
==
2
:
v
=
np
.
insert
(
v
,
self
.
_dim
,
0
,
axis
=
1
)
da
=
np
.
insert
(
da
,
self
.
_dim
,
0
,
axis
=
1
)
us
=
np
.
empty
((
self
.
n_nodes
(),
self
.
dimension
()))
self
.
_lib
.
fluid_problem_u_solid
(
self
.
_fp
,
_np2c
(
us
))
data
=
[
(
"pressure"
,
self
.
solution
()[:,[
self
.
_dim
]]),
(
"velocity"
,
v
),
(
"porosity"
,
self
.
porosity
()),
(
"old_porosity"
,
self
.
old_porosity
()),
(
"us"
,
us
),
(
"grad"
,
da
),
(
"parent_node_id"
,
self
.
parent_nodes
())
]
...
...
@@ -490,7 +495,7 @@ class FluidProblem :
contact -- List of particles contact resultant forces
"""
self
.
n_particles
=
mass
.
shape
[
0
]
self
.
_lib
.
fluid_problem_set_particles
(
self
.
_fp
,
c_int
(
mass
.
shape
[
0
]),
_np2c
(
mass
),
_np2c
(
volume
),
_np2c
(
position
),
_np2c
(
velocity
),
_np2c
(
contact
))
self
.
_lib
.
fluid_problem_set_particles
(
self
.
_fp
,
c_int
(
mass
.
shape
[
0
]),
_np2c
(
mass
),
_np2c
(
volume
),
_np2c
(
position
),
_np2c
(
velocity
),
_np2c
(
contact
))
def
move_particles
(
self
,
position
,
velocity
,
contact
):
"""Set location of the grains in the mesh and compute the porosity in each cell.
...
...
python/time_integration.py
View file @
8591a6ba
...
...
@@ -127,7 +127,7 @@ def predictor_corrector_iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e
# Give to the fluid the solid information
fluid
.
move_particles
(
particles
.
position
(),
particles
.
velocity
(),
particles
.
contact_forces
())
def
iterate
(
fluid
,
particles
,
dt
,
min_nsub
=
1
,
contact_tol
=
1e-8
,
external_particles_forces
=
None
,
fixed_grains
=
False
,
after_sub_iter
=
None
,
max_nsub
=
None
,
check_residual_norm
=-
1
)
:
def
iterate
(
fluid
,
particles
,
dt
,
min_nsub
=
1
,
contact_tol
=
1e-8
,
external_particles_forces
=
None
,
fixed_grains
=
False
,
after_sub_iter
=
None
,
max_nsub
=
None
,
check_residual_norm
=-
1
,
use_predictor_corrector
=
True
)
:
"""Migflow solver: the solver type depends on the fluid and particles given as arguments.
Keyword arguments:
...
...
@@ -161,10 +161,20 @@ def iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e-8, external_particl
if
fluid
is
not
None
and
particles
is
not
None
and
particles
.
r
()
is
not
None
and
particles
.
r
().
size
!=
0
:
if
fixed_grains
:
fluid
.
implicit_euler
(
dt
,
check_residual_norm
)
if
fluid
.
n_fluids
()
==
2
:
fluid
.
advance_concentration
(
dt
)
fluid
.
move_particles
(
particles
.
position
(),
particles
.
velocity
(),
-
fluid
.
compute_node_force
(
dt
))
else
:
predictor_corrector_iterate
(
fluid
,
particles
,
dt
,
min_nsub
,
contact_tol
,
external_particles_forces
,
after_sub_iter
=
after_sub_iter
,
max_nsub
=
max_nsub
,
check_residual_norm
=
check_residual_norm
)
if
use_predictor_corrector
:
predictor_corrector_iterate
(
fluid
,
particles
,
dt
,
min_nsub
,
contact_tol
,
external_particles_forces
,
after_sub_iter
=
after_sub_iter
,
max_nsub
=
max_nsub
,
check_residual_norm
=
check_residual_norm
)
else
:
fluid
.
implicit_euler
(
dt
,
check_residual_norm
)
f0
=
fluid
.
compute_node_force
(
dt
)
+
external_particles_forces
_advance_particles
(
particles
,
f0
,
dt
,
min_nsub
,
contact_tol
,
max_nsub
=
max_nsub
)
if
fluid
.
n_fluids
()
==
2
:
fluid
.
advance_concentration
(
dt
)
fluid
.
move_particles
(
particles
.
position
(),
particles
.
velocity
(),
particles
.
contact_forces
())
class
FreeSurface
():
...
...
testcases/depot-2d/depot.py
View file @
8591a6ba
...
...
@@ -75,7 +75,7 @@ rho = 1000 # fluid density
nu
=
1e-6
# kinematic viscosity
# Numerical parameters
outf
=
5
# number of iterations between output files
outf
=
100
# number of iterations between output files
dt
=
1e-3
# time step
tEnd
=
100
# final time
...
...
@@ -102,7 +102,7 @@ ii = 0
#
# FLUID PROBLEM
#
fluid
=
fluid
.
FluidProblem
(
2
,
g
,[
nu
*
rho
],[
rho
])
fluid
=
fluid
.
FluidProblem
(
2
,
g
,[
nu
*
rho
],[
rho
]
,
drag_in_stab
=
0
,
solver
=
"petsc"
,
solver_options
=
"-pc_type lu"
,
usolid
=
True
)
# Set the mesh geometry for the fluid computation
fluid
.
load_msh
(
"mesh.msh"
)
fluid
.
set_wall_boundary
(
"Bottom"
)
...
...
@@ -118,7 +118,7 @@ tic = time.time()
# COMPUTATION LOOP
#
while
t
<
tEnd
:
time_integration
.
iterate
(
fluid
,
p
,
dt
,
min_nsub
=
5
,
external_particles_forces
=
g
*
p
.
mass
())
time_integration
.
iterate
(
fluid
,
p
,
dt
,
min_nsub
=
5
,
external_particles_forces
=
g
*
p
.
mass
()
,
use_predictor_corrector
=
False
)
t
+=
dt
# Output files writting
...
...
validation/darcy/darcy.py
View file @
8591a6ba
...
...
@@ -108,9 +108,8 @@ class Darcy(unittest.TestCase) :
fluid
=
mbfluid
.
FluidProblem
(
2
,
g
,
mu
,
rho
,
drag_in_stab
=
True
)
fluid
.
load_msh
(
"mesh.msh"
)
U
=
0.00001
# averaged fluid velocity
V
=
U
/
(
1
-
compacity
)
# fluid velocity
fluid
.
solution
()[:,
0
]
=
U
fluid
.
set_open_boundary
(
"Left"
,
velocity
=
[
V
,
0
])
fluid
.
set_open_boundary
(
"Left"
,
velocity
=
[
U
,
0
])
fluid
.
set_open_boundary
(
"Right"
,
pressure
=
0
)
fluid
.
set_wall_boundary
(
"Bottom"
,
velocity
=
[
0
,
0
])
fluid
.
set_wall_boundary
(
"Top"
,
velocity
=
[
0
,
0
])
...
...
Write
Preview
Supports
Markdown
0%
Try again
or
attach a new file
.
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Cancel
Please
register
or
sign in
to comment