Commit 92f33f44 authored by Matthieu Constant's avatar Matthieu Constant
Browse files

new automatic testcases

parent 5a7790ca
Pipeline #3733 failed with stage
in 14 seconds
# 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 mbfluid
from marblesbag import scontact2
import numpy as np
import os
import subprocess
import time
import shutil
import random
import unittest
#Physical parameters for the drops are the ones presented by Metzger et al. (2007) "Falling clouds of particles in viscous fluids"
class Poiseuille(unittest.TestCase) :
def runTest(self) :
dir_path = os.path.dirname(os.path.realpath(__file__))
os.chdir(dir_path)
# 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)
subprocess.call(["gmsh", "-2", "mesh.geo","-clscale","2"])
t = 0
ii = 0
#physical parameters
g = 0 # 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 = 0.1 # time step
alpha = 1e-6 # stabilization coefficient
epsilon = alpha*lcmin**2 /nu # stabilization parametre
print('epsilon',epsilon)
shutil.copy("mesh.msh", outputdir +"/mesh.msh")
vUp = 1
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,vUp),
("Top",1,1,0),
("PtFix",3,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 = mbfluid.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()
p = [ 4.40910528e+03, -2.08531786e+04, 4.17388709e+04, -4.57353622e+04, 2.95421532e+04, -1.12048001e+04 ,2.30784236e+03 ,-2.14192853e+02, 1.42357618e+01, -3.69376852e+00, -5.13246798e-04]
E = 0
N = 0
for i in range(len(x)):
if abs(x[i,0])<0.01:
E += abs(s[i,0]-np.polyval(p,x[i,1]))
N += 1
self.assertLess(E/N,5*vUp/100, "error is too large in Cavity")
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") = {2};
import unittest
all_tests = [
"poiseuille.poiseuille",
"poiseuille.poiseuille-2fluids"]
"poiseuille.poiseuille-2fluids",
"cavity.cavity"]
suite = unittest.TestSuite()
for t in all_tests :
suite.addTest(unittest.defaultTestLoader.loadTestsFromName(t))
......
......@@ -102,10 +102,12 @@ class Poiseuille(unittest.TestCase) :
#set initial_condition
x = fluid.coordinates()
s = fluid.solution()
a = (x[:,0]+5)/10
s[:,2] = x[:,1]
s[:,5] = 1-s[:,2]
mu = nu0*rho
vI = 1/(120*mu)
fluid.export_vtk(outputdir,0,0)
ii = 0
......@@ -124,6 +126,8 @@ class Poiseuille(unittest.TestCase) :
s = fluid.solution()
x = fluid.coordinates()
vel = (s[:,0]+s[:,3]-1/(20*nu0*rho)*x[:,1]*(1-x[:, 1]))**2
vS = np.sum((1/(20*nu0*rho)*x[:,1]*(1-x[:, 1]))**2)
print('Error', (vel.sum())**.5)
self.assertLess(vel.sum()**.5,2.2e-3, "error is too large in Poiseuille 2 fluids")
self.assertLess(vel.sum()**.5,(vS**0.5)/100, "error is too large in Poiseuille")
......@@ -74,10 +74,8 @@ class Poiseuille(unittest.TestCase) :
#set initial_condition
x = fluid.coordinates()
s = fluid.solution()
a = (x[:,0]+5)/10
s[:,2] = 1.
fluid.export_vtk(outputdir,0,0)
ii = 0
......@@ -96,6 +94,7 @@ class Poiseuille(unittest.TestCase) :
s = fluid.solution()
x = fluid.coordinates()
vel = (s[:,0]-1/(20*nu*rho)*x[:,1]*(1-x[:, 1]))**2
vS = np.sum((1/(20*nu*rho)*x[:,1]*(1-x[:, 1]))**2)
print('Error', (vel.sum())**.5)
self.assertLess(vel.sum()**.5,1.5e-3, "error is too large in Poiseuille")
self.assertLess(vel.sum()**.5,(vS**0.5)/100, "error is too large in Poiseuille")
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