Commit c6bbcbb9 authored by Michel Henry's avatar Michel Henry
Browse files

test

parent dcb72f70
Pipeline #8311 passed with stages
in 3 minutes and 9 seconds
......@@ -2,22 +2,22 @@
* MigFlow - Copyright (C) <2010-2018>
* <Universite catholique de Louvain (UCL), Belgium
* Universite de Montpellier, France>
*
*
* List of the contributors to the development of MigFlow: see AUTHORS file.
* Description and complete License: see LICENSE file.
*
* This program (MigFlow) is free software:
* you can redistribute it and/or modify it under the terms of the GNU Lesser General
*
* This program (MigFlow) 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,
* along with this program (see COPYING and COPYING.LESSER files). If not,
* see <http://www.gnu.org/licenses/>.
*/
......@@ -61,9 +61,9 @@ inline static double mesh_dxidx(const Mesh *mesh, int iel, double dxidx[D][D]) {
inline static double element_volume_from_detj(double detj) {
#if D == 2
return detj/2;
#else
#else
return detj/6;
#endif
#endif
}
inline static double node_volume_from_detj(double detj) {
......@@ -347,7 +347,7 @@ static void particle_force_f(FluidProblem *problem, double *f, double *dfdu, dou
massreduced = (rhop-problem->rho[0])*problem->particle_volume[ip];
}
double gamma;
if(problem->n_fluids == 2){
double gamma0 = particle_drag_coeff(Due,problem->mu[0],problem->rho[0],vol,c);
double gamma1 = particle_drag_coeff(Due,problem->mu[1],problem->rho[1],vol,c);
......@@ -396,7 +396,7 @@ static void fluid_problem_f(const FluidProblem *problem, const double *sol, doub
double rhoreduced = (problem->reduced_gravity ? (rho-problem->rho[0]) : rho);
double drag = problem->volume_drag + problem->quadratic_drag*nold;
for (int i = 0; i < D; ++i) {
f0[U+i] =
f0[U+i] =
+ c*dp[i]
- g[i]*rhoreduced*c
- bf[i]*c
......@@ -422,13 +422,13 @@ static void fluid_problem_f(const FluidProblem *problem, const double *sol, doub
f1[(U+i)*D+j] = mu*(du[i][j]-u[i]*dc[j]/c+ du[j][i]-u[j]*dc[i]/c);
f10[((U+i)*D+j)*n_fields+U+j] += -mu*dc[i]/c;
f10[((U+i)*D+j)*n_fields+U+i] += -mu*dc[j]/c;
f11[((U+i)*D+j)*n_fields*D+(U+j)*D+i] += mu;
f11[((U+i)*D+j)*n_fields*D+(U+j)*D+i] += mu;
f11[((U+i)*D+j)*n_fields*D+(U+i)*D+j] += mu;
// SUPG
double supg = uold[j]*taup;
f1[(U+i)*D+j] += supg*f0[U+i];
f10[((U+i)*D+j)*n_fields+U+i] += supg*f00[(U+i)*n_fields+U+i];
f11[((U+i)*D+j)*(n_fields*D)+P*D+i] += supg*f01[(U+i)*n_fields*D+P*D+i];
f11[((U+i)*D+j)*(n_fields*D)+P*D+i] += supg*f01[(U+i)*n_fields*D+P*D+i];
for (int k = 0; k < D; ++k)
f11[((U+i)*D+j)*n_fields*D+(U+i)*D+k] += supg*f01[(U+i)*n_fields*D+(U+i)*D+k];
}
......@@ -610,7 +610,7 @@ double fluid_problem_a_integ_volume(FluidProblem *problem)
double dxidx[D][D];
const double detj = mesh_dxidx(mesh,iel,dxidx);
for (int i = 0; i < N_N; ++i)
for (int j = 0; j< N_N; ++j)
for (int j = 0; j< N_N; ++j)
I_a += detj * problem->porosity[el[i]]*problem->concentration[iel*N_N+j]*mass_matrix[i][j];
}
return I_a;
......@@ -704,7 +704,7 @@ double fluid_problem_a_integ_bnd(FluidProblem *problem, double dt)
double n[D],detbnd;
get_normal_and_det(x,n,&detbnd);
grad_shape_functions(dxidx, dphi);
for (int iqp = 0; iqp < N_LQP; ++iqp) {
double *dataqp = &data[n_value*(N_LQP*iinterface+iqp)];
double phi[DIMENSION];
......@@ -1120,7 +1120,7 @@ void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, const dou
const Mesh *mesh = problem->mesh;
int n_fields = fluid_problem_n_fields(problem);
size_t local_size = N_SF*n_fields;
size_t N_E = mesh->n_elements;
for (size_t i = 0; i < N_E*local_size; ++i)
......@@ -1142,6 +1142,7 @@ void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, const dou
printf("Boundary tag \"%s\" not found.\n", bnd->tag);
exit(1);
}
// To change for periodic conditions
for (int i = 0; i < problem->boundaries[iphys].n_nodes; ++i){
forced_row[problem->boundaries[iphys].nodes[i]*n_fields+bnd->row] = bnd->field;
}
......@@ -1163,14 +1164,18 @@ void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, const dou
const int *el = &mesh->elements[iel*N_N];
char forced_field = forced_row[el[inode]*n_fields+ifield];
if (forced_field == -1) continue;
int i = inode*n_fields + ifield;
int i = inode*n_fields + ifield; // HERE : To periodic
// int i = parent_tag[inode]*n_fields + ifield;
for (int jnode = 0; jnode< N_SF; ++jnode){
for (int jfield= 0; jfield< n_fields; ++jfield){
// int j = parent_tag[jnode]*n_fields+jfield
local_matrix[(inode*n_fields+ifield)*local_size+(jnode*n_fields+jfield)] =
(inode==jnode && jfield == forced_field) ? 1. : 0.;
// local_matrix[i*local_size+j] = (inode==jnode && jfield == forced_field) ? 1. : 0.;
}
}
local_vector[inode+ifield*N_SF] = 0;
// local_vector[parent_tag[inode]+ifield*N_SF] = 0;
}
}
}
......@@ -1464,7 +1469,7 @@ void fluid_problem_adapt_gen_mesh(FluidProblem *problem, double lcmax, double lc
double lc = lcmin /fmin(1,fmax(ngrad/(gradmax),gradmin/(gradmax)));
if (problem->n_fluids == 1){
double ngradP = 0.;
double ngradP = 0.;
for (int j=0; j<DIMENSION; ++j){
double gradP = 0;
for (int k=0; k<N_SF; ++k){
......@@ -1476,7 +1481,7 @@ void fluid_problem_adapt_gen_mesh(FluidProblem *problem, double lcmax, double lc
double lcp = lcmin /fmin(1,fmax(ngradP2/(gradPmax),gradPmin/(gradPmax)));
lc = fmin(lcp,lc);
}
if(problem->n_fluids == 2){
double ngradA = 0;
for (int k = 0; k < D; ++k){
......@@ -1646,7 +1651,7 @@ void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin,
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 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(problem->mesh_tree,new_mesh->n_nodes, newx, new_eid, new_xi);
......@@ -1730,10 +1735,10 @@ void fluid_problem_after_import(FluidProblem *problem)
void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *state, double *contact, long int *elid, int init, int reload)
{
struct timespec tic, toc;
int n_fields = fluid_problem_n_fields(problem);
if(problem->n_particles != n) {
problem->n_particles = n;
free(problem->particle_uvw);
......@@ -1967,7 +1972,7 @@ void weak_boundary_values(const Mesh *mesh, const MeshBoundary *bnd, const WeakB
for (int iphi = 0; iphi < D; ++iphi) {
xx += mesh->x[el[cl[iphi]]*3+id]*phil[iphi];
}
x[(i*N_LQP+iqp)*D+id] = xx;
x[(i*N_LQP+iqp)*D+id] = xx;
}
}
}
......@@ -1977,6 +1982,6 @@ void weak_boundary_values(const Mesh *mesh, const MeshBoundary *bnd, const WeakB
}
void fluid_problem_reset_boundary_force(FluidProblem *problem) {
for (int i=0;i<D*problem->mesh->n_boundaries;i++)
for (int i=0;i<D*problem->mesh->n_boundaries;i++)
problem->boundary_force[i] = 0;
}
L = 5;
H = 1;
y = 0;
lc = 1;//.05;
lc = 0.05;
Point(1) = {0,0,0};
Point(2) = {1,0,0};
Point(3) = {1,1,0};
Point(4) = {0,1,0};
Point(1) = {0,0,0,lc};
Point(2) = {1,0,0,lc};
Point(3) = {1,1,0,lc};
Point(4) = {0,1,0,lc};
Line(1) = {1,2};
Line(2) = {2,3};
Line(3) = {4,3};
Line(4) = {1,4};
Line Loop(1) = {1,2,-3,-4};
Plane Surface(1) = {1};
Periodic Curve {1} = {3};
Periodic Curve {1,2} = {3,4};
Physical Line("Left") = {4};
Physical Line("Right") = {2};
......
......@@ -48,7 +48,7 @@ g = np.array([0,0]) # gravity
rho = 1000 # fluid density
nu = 1e-3 # kinematic viscosity
mu = nu*rho # dynamic viscosity
tEnd = 100000 # final time
tEnd = 10000 # final time
#numerical parameters
lcmin = .1 # mesh size
......@@ -65,7 +65,7 @@ outf = 1 # number of iterations between ou
fluid = mbfluid.FluidProblem(2,g,nu*rho,rho)
fluid.load_msh("mesh.msh")
fluid.set_open_boundary("Left",velocity=[lambda x : 1/(20*mu)*x[:,1]*(1-x[:, 1]),0])
fluid.set_wall_boundary("Bottom",velocity=[0,0])
fluid.set_wall_boundary("Bottom",velocity=[0,0], pressure=0)
fluid.set_wall_boundary("Top",velocity=[0,0], pressure=0)
fluid.set_open_boundary("Right",velocity=[lambda x : 1/(20*mu)*x[:,1]*(1-x[:, 1]),0])
......@@ -79,7 +79,7 @@ fluid.export_vtk(outputdir,0,0)
ii = 0
tic = time.time()
while ii < 100 :
while t < tEnd :
#Fluid solver
time_integration.iterate(fluid,None,dt)
t += dt
......
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