Commit 4595d5ca authored by Matthieu Constant's avatar Matthieu Constant
Browse files

gravite + termes de bord + cas test gravite

parent 7f50552a
......@@ -252,24 +252,48 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, d
const int nodes[2] = {el[i0],el[i1]};
const double *x[2] = {&mesh->x[nodes[0]*3], &mesh->x[nodes[1]*3]};
const double dx[2] = {x[1][0]-x[0][0], x[1][1]-x[0][1]};
double dxdxi[D][D], dxidx[D][D], dphi[N_N][D];
for (int i = 0; i < D; ++i)
for (int j = 0; j < D; ++j)
dxdxi[i][j] = mesh->x[el[j+1]*3+i] - mesh->x[el[0]*3+i];
const double detj = invDD(dxdxi, dxidx);
grad_shape_functions(dxidx, dphi);
double dp[D]={0};
const double l = sqrt(dx[0]*dx[0] + dx[1]*dx[1]);
const double n[2] = {-dx[1]/l,dx[0]/l};
double *local_matrix = &all_local_matrix[local_size*local_size*iel];
double *local_vector = &all_local_vector[local_size*iel];
const int U = 0;
const int Q = 2;
const int P = 3;
for (int i = 0; i< N_SF; ++i) {
for (int j = 0; j < D; ++j) {
dp[j] += dphi[i][j]*solution[el[i]*n_fields+3];
}
}
double p0 = solution[el[i0]*n_fields+P];
double p1 = solution[el[i1]*n_fields+P];
/*for (int k = 0; k < D; ++k) {
for (int k = 0; k < D; ++k) {
if (!forced) {
local_matrix[(i0*n_fields+U+k)*local_size + (i0*n_fields+P)] += l/3*n[k];
local_matrix[(i0*n_fields+U+k)*local_size + (i1*n_fields+P)] += l/6*n[k];
local_matrix[(i1*n_fields+U+k)*local_size + (i0*n_fields+P)] += l/6*n[k];
local_matrix[(i1*n_fields+U+k)*local_size + (i1*n_fields+P)] += l/3*n[k];
local_vector[i0*n_fields+U+k] += l*(p0/3+p1/6)*n[k];
local_vector[i1*n_fields+U+k] += l*(p0/6+p1/3)*n[k];
local_vector[i0+(U+k)*N_SF] += l*(p0/3+p1/6)*n[k];
local_vector[i1+(U+k)*N_SF] += l*(p0/6+p1/3)*n[k];
local_matrix[(i0*n_fields+Q)*local_size + (i0*n_fields+P)] += dphi[i0][k]*n[k]*problem->epsilon/2.*l;
local_matrix[(i0*n_fields+Q)*local_size + (i1*n_fields+P)] += dphi[i1][k]*n[k]*problem->epsilon/2.*l;
local_matrix[(i1*n_fields+Q)*local_size + (i0*n_fields+P)] += dphi[i0][k]*n[k]*problem->epsilon/2.*l;
local_matrix[(i1*n_fields+Q)*local_size + (i1*n_fields+P)] += dphi[i1][k]*n[k]*problem->epsilon/2.*l;
local_vector[i0+Q*N_SF] += dp[k]*n[k]*problem->epsilon/2.*l;
local_vector[i1+Q*N_SF] += dp[k]*n[k]*problem->epsilon/2.*l;
}
}*/
}
if (forced) {
for (int ifluid = 0; ifluid < problem->n_fluids; ++ifluid) {
const int Q = ifluid*(D+1)+D;
......@@ -298,8 +322,8 @@ static void fluid_problem_f(FluidProblem *problem, double *f0, double *f1, doubl
double epsilonp = problem->epsilon;
double epsilonq = problem->epsilon/10;
f0[P] = 1-c;
f1[P*D+0] = epsilonp*dp[0];
f1[P*D+1] = epsilonp*dp[1];
f1[P*D+0] = 0;//epsilonp*dp[0];
f1[P*D+1] = 0;//epsilonp*dp[1];
for (int ifluid=0; ifluid < problem->n_fluids; ++ifluid) {
int Q = ifluid*(D+1)+D;
int U = ifluid*(D+1);
......@@ -327,11 +351,11 @@ static void fluid_problem_f(FluidProblem *problem, double *f0, double *f1, doubl
f0[Q] = divu+(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];
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);
for (int j = 0; j < D; ++j) {
f1[(U+j)*D+i] = mu*(tau[i][j]+tau[j][i]) + (i==j ? -q*p : 0);
}
f1[Q*D+i] = epsilonq*dq[i];
f1[Q*D+i] = epsilonq*dq[i]+epsilonp*q*dp[i];
}
}
}
......@@ -465,6 +489,7 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
}
}
}
}
for (int inode = 0; inode < N_SF; ++inode) {
for(int ifield = 0; ifield < n_fields; ++ifield) {
......
L = .5;
H = 1;
y = 0;
lc = .05;
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};
Line Loop(1) = {1:4};
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") = {1};
# 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 = 10 # time step
alpha = 1e-5 # 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",1,1,0.),("Bottom",1,1,0.),("Right",1,1,0.),("Left",1,1,0.),
("Top",0,0,0.),("Bottom",0,0,0.),("Left",0,0,0),("Right",0,0,0),
("Top",2,3,0)]#,("Bottom",2,3,-rho*g)]#,("Left",2,3,lambda x : -rho*g*(1-x[:,1])),("Right",2,3,lambda x : -rho*g*(1-x[:,1]))]
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)
\ No newline at end of file
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