Commit 5d5c6f69 authored by Michel Henry's avatar Michel Henry
Browse files

wip

parent 6e39905b
Pipeline #9978 failed with stages
in 2 minutes and 41 seconds
......@@ -414,22 +414,19 @@ static void fluid_problem_f(const FluidProblem *problem, const double *sol, doub
double *g = problem->g;
double rhoreduced = (problem->reduced_gravity ? (rho-problem->rho[0]) : rho);
double drag = problem->volume_drag + problem->quadratic_drag*nold;
// f0[P] = p;
// f00[P*n_fields+P] = 1;
for (int i = 0; i < D; ++i) {
f0[U+i] =
+ c*dp[i]
- g[i]*rhoreduced*c
- bf[i]*c
+ drag*u[i];
+ 5.3e5*u[i];
// + 5.3e5*u[i];
f00[(U+i)*n_fields+U+i] += drag;//5.3e5;
f01[(U+i)*n_fields+P][i] += c;
if(problem->temporal) {
f0[U+i] += rho*(u[i]-uold[i])/dt;// - rho;
f0[U+i] += rho*(u[i]-uold[i])/dt;
f00[(U+i)*n_fields+U+i] += rho/dt;
}
// #if 0
if (problem->advection) {
// advection :
// dj(uj ui /c) = uj/c dj(ui) + ui/c dj(uj) - uj ui /(c*c) dj(c)
......@@ -472,7 +469,6 @@ static void fluid_problem_f(const FluidProblem *problem, const double *sol, doub
for (int j = 0; j < D; ++j) {
f11[P*n_fields+U+i][i][j] = pspg*f01[(U+i)*n_fields+(U+i)][j];
}
// #endif
}
}
......@@ -1110,7 +1106,8 @@ void fluid_problem_apply_boundary_conditions(FluidProblem *problem)
double *x = malloc(ndof_closure*sizeof(double)*D);
double *v = malloc(ndof_closure*sizeof(double));
for(int id=0; id < ndof_closure; ++id){
int closure = (*field->closure + ndof_closure*interface[1])[id];
// int closure = field->closure + ndof_closure*interface[1])[id];
int closure = field->closure[interface[1]*ndof_closure+id];
double xibnd[D];
for(int d=0; d < D; ++d)
xibnd[d] = field->param_coord[closure][d];
......@@ -1133,7 +1130,7 @@ void fluid_problem_apply_boundary_conditions(FluidProblem *problem)
}
int *dofmap = &(map[imap]);
for(int id=0; id < ndof_closure; ++id){
int closure = (*field->closure + ndof_closure*interface[1])[id];
int closure = field->closure[interface[1]*ndof_closure+id];
problem->solution[dofmap[closure]] = v[id];
}
free(x);
......
......@@ -4,13 +4,8 @@
#include <string.h>
static int CLOSURE_TRIANGLE_P1[6][2] = {
{0,1},{1,2},{2,0},
{1,0},{2,1},{0,2}};
static int CLOSURE_TRIANGLE_P2[6][3] = {
{0,1,3},{1,2,4},{2,0,5},
{1,0,3},{2,1,4},{0,2,5}};
static int CLOSURE_TRIANGLE_P1[12] = {0,1,1,2,2,0,1,0,2,1,0,2};
static int CLOSURE_TRIANGLE_P2[18] = {0,1,3,2,2,4,2,0,5,1,0,3,2,1,4,0,2,5};
static double XI_TRIANGLE_P1[][D] = {
{0.0,0.0},{1.0,0.0},{0.0,1.0}};
......@@ -314,8 +309,7 @@ void fe_fields_add_to_local_vector(const FEFields *fields, const double *f0, con
v += dsf[cc][id]*f1[ifield][id]*jw;
}
}
local_vector[cc] += v;
cc++;
local_vector[cc++] += v;
}
}
}
......
......@@ -9,7 +9,7 @@ struct FEElementStruct{
int nlocal;
int n[DIMENSION+1];
int ndof_closure;
int (*closure)[];
int *closure;
double (*param_coord)[D];
void (*f)(const double *xi, double *f);
void (*df)(const double *xi, const double dxidx[D][D], double f[][D]);
......@@ -71,56 +71,23 @@ static void fe_fields_add_to_local_matrix(const FEFields *fields, const double *
ic++;
}
}
// int c = 0;
// int jc = 0;
// int n_fields = fields->n;
// for (int jfield = 0; jfield < fields->n; ++jfield) {
// for (int jphi = 0; jphi < fields->element[jfield]->nlocal; ++jphi){
// int ic = 0;
// for (int ifield = 0; ifield < fields->n; ++ifield){
// for (int iphi = 0; iphi < fields->element[ifield]->nlocal; ++iphi){
// double m = 0;
// if (f00) {
// m += jw*sf[ic]*sf[jc]*f00[ifield*fields->n+jfield];
// }
// if (f10) {
// for (int id = 0; id < D; ++id) {
// m += jw*dsf[ic][id]*sf[jc]*f10[ifield*fields->n+jfield][id];
// }
// }
// if (f11) {
// for (int id = 0; id < D; ++id) {
// for (int jd = 0; jd < D; ++jd) {
// m += jw*dsf[ic][id]*dsf[jc][jd]*f11[ifield*fields->n+jfield][id][jd];
// }
// }
// }
// if (f01) {
// for (int jd = 0; jd < D; ++jd) {
// m += jw*sf[ic]*dsf[jc][jd]*f01[ifield*fields->n+jfield][jd];
// }
// }
// local_matrix[c++] += m;
// ic ++;
// }
// }
// jc ++;
// }
// }
}
void fe_fields_free(FEFields *fields);
#if D==2
// #define N_LQP 2
// #define N_QP 3
#define N_LQP 4
#define N_QP 13
#define N_LQP 2
static double LQP[][1] = {{0.21132486540518713}, {0.7886751345948129}};
static double LQW[] = {0.5,0.5};
// static double LQP[][1] = {{0.21132486540518713}, {0.7886751345948129}};
// static double LQW[] = {0.5,0.5};
static double LQW[] = {0.173927422568726842,0.326072577431273214,0.326072577431273214,0.173927422568726842};
static double LQP[][1] = {{0.0694318442029737137},{0.330009478207571871},{0.669990521792428129},{0.930568155797026231}};
// static double QP[][2] = {{1./6,1./6},{1./6,2./3},{2./3,1./6}};
// static double QW[] = {1./6,1./6,1./6};
static double QW[] = {-0.0747850222338409948,0.087807628716603997,0.087807628716603997,0.087807628716603997,0.0266736178044190003,0.0266736178044190003,0.0266736178044190003,0.0385568804451279987,0.0385568804451279987,0.0385568804451279987,0.0385568804451279987,0.0385568804451279987,0.0385568804451279987};
static double QP[][2] = {{0.333333333333332982,0.333333333333332982},{0.260345966079039981,0.260345966079039981},{0.260345966079039981,0.479308067841919982},{0.479308067841919982,0.260345966079039981},{0.0651301029022160055,0.0651301029022160055},{0.0651301029022160055,0.869739794195568017},{0.869739794195568017,0.0651301029022160055},{0.312865496004874011,0.638444188569810001},{0.638444188569810001,0.0486903154253160025},{0.0486903154253160025,0.312865496004874011},{0.312865496004874011,0.0486903154253160025},{0.638444188569810001,0.312865496004874011},{0.0486903154253160025,0.638444188569810001}};
......
import numpy as np
from ._tools import gmsh
# import gmsh
# gmsh.initialize()
#static p1 arrays, could be generated automatically but hard-coded to match gmsh order
......
......@@ -356,21 +356,18 @@ class FluidProblem :
t -- Computation time
stab -- If True exports the stabilisation parametres in the output files
"""
v = self.solution()[:,:self._dim]
da = self.concentration_dg_grad()
cell_data = []
if self._dim == 2 :
v = np.insert(v,self._dim,0,axis=1)
da = np.insert(da,self._dim,0,axis=1)
data = [
("pressure",self.solution()[:,[self._dim]]),
("velocity",v),
("pressure",self.pressure().reshape(-1,1)),
("porosity",self.porosity()),
("old_porosity",self.old_porosity()),
("grad",da),
("parent_node_id", self.parent_nodes())
]
if self._n_fluids == 2 :
data.append(("grad",da))
cell_data.append(("curvature",self.curvature()))
cell_data.append(("concentration",self.concentration_dg()))
cell_data.append(("stab",self._get_matrix("stab_param",self.n_elements(),1)))
......@@ -380,7 +377,29 @@ class FluidProblem :
offsets = np.cumsum(np.repeat([self._dim+1],connectivity.shape[0])).astype(np.int32)
VTK.write(output_dir+"/fluid",it,t,(types,connectivity,offsets),self.coordinates(),data,field_data,cell_data)
def read_vtk(self, dirname, i): # to update for p2p1
v = self.velocity()
if self._dim == 2 :
v = np.insert(v,self._dim,0,axis=1)
data = [("velocity",v)]
if np.max(self.order_fields) == 1:
connectivity = self.elements()
types = np.repeat([5 if self._dim == 2 else 10],connectivity.shape[0])
offsets = np.cumsum(np.repeat([self._dim+1],connectivity.shape[0])).astype(np.int32)
VTK.write(output_dir+"/velocity_p1",it,t,(types,connectivity,offsets),self.coordinates(),data,None,None)
elif np.max(self.order_fields) == 2:
connectivity = self.set_numbering([2]).reshape(self.n_elements(),-1)
types = np.repeat([22 if self._dim == 2 else 24],connectivity.shape[0])
offsets = np.cumsum(np.repeat([connectivity.shape[1]],connectivity.shape[0])).astype(np.int32)
xp2 = self.coordinates_p2()
VTK.write(output_dir+"/velocity_p2",it,t,(types,connectivity,offsets),xp2,data,None,None)
else:
raise ValueError("Order unavailable !")
def read_vtk(self, dirname, i):
"""Reads output file to reload computation.
Keyword arguments:
......@@ -403,9 +422,14 @@ class FluidProblem :
_np2c(bnd_tags,np.int32),c_int(len(cbnd_names)),cbnd_names,
_np2c(data["parent_node_id"],np.int32) if "parent_node_id" in data else None))
self._lib.fluid_problem_set_mesh(self._fp, _mesh)
ndof = np.max(self.set_numbering(self.order_fields))+1
id_v = np.r_[np.repeat(np.arange(self._dim)[None,:],self.n_nodes(),axis=0).reshape(-1)
+ self.n_fields()*np.repeat(np.arange(self.n_nodes()), self._dim),
np.arange(self.n_fields()*self.n_nodes(), ndof)]
id_p = np.arange(self._dim, self.n_fields()*self.n_nodes(), self.n_fields())
sol = self.solution()
sol[:,:self._dim] = data["velocity"][:,:self._dim]
sol[:,[self._dim]] = data["pressure"]
# sol[id_v] = data["velocity"][:,:self._dim].reshape([-1])
sol[id_p] = data["pressure"].reshape([-1])
if self._n_fluids == 2 :
self.concentration_dg()[:] = cdata["concentration"]
self.porosity()[:] = data["porosity"]
......@@ -456,11 +480,9 @@ class FluidProblem :
self._lib.fluid_problem_porosity_local_map(self._fp, _np2c(idx_porosity, dtype=np.int32))
self._lib.fluid_problem_concentration_local_map(self._fp, _np2c(idx_concentration, dtype=np.int32))
# print(time()-timer)
# idx = np.zeros((self.n_elements()*self.local_size()), dtype=np.int32)
# self.local_map(idx)
idx = idx.reshape([self.n_elements(), -1])
periodic_map = self._get_matrix("periodic_mapping", self.n_nodes(), 1, c_int).flatten()
solution_old = self.solution_flat().copy()
solution_old = self.solution().copy()
if self.sys is None :
constraints = []
if self.mean_pressure is not None:
......@@ -475,19 +497,14 @@ class FluidProblem :
localm = np.ndarray([sys.localsize**2*self.n_elements()])
self._lib.fluid_problem_assemble_system(self._fp,_np2c(solution_old),c_double(dt),_np2c(localv),_np2c(localm))
# print(time()-timer)
# np.set_printoptions(precision=1,threshold=1e-8)
# mass_matrix = localm[:sys.localsize**2].reshape(sys.localsize, -1)
# print(360*mass_matrix[:6,:6])
constraint_value = []
if self.mean_pressure is not None:
constraint_value.append((self.node_volume().reshape([-1]), -self.mean_pressure*np.sum(self.node_volume())))
rhs = sys.local_to_global(localv,localm,self.solution_flat(),constraint_value)
rhs = sys.local_to_global(localv,localm,self.solution(),constraint_value)
# self.solution()[:] -= sys.solve(rhs).reshape([-1, self.n_fields()])[periodic_map]
# self.solution()[:] -= sys.solve(rhs).reshape([-1, self.n_fields()])
self.solution_flat()[:] -= sys.solve(rhs)
# exit(0)
self.solution()[:] -= sys.solve(rhs)
# print(time()-timer)
if check_residual_norm > 0 and False :
if check_residual_norm > 0:
self._lib.fluid_problem_assemble_system(self._fp,_np2c(solution_old),c_double(dt),_np2c(localv),_np2c(localm))
rhs = sys.local_to_global(localv,localm,self.solution().reshape([-1]),constraint_value)
norm = np.linalg.norm(rhs)
......@@ -501,7 +518,7 @@ class FluidProblem :
"""Computes the CFL number divided by the time step
"""
nv = np.linalg.norm(self.solution()[:,:self._dim],axis=1)
nv = np.linalg.norm(self.velocity(),axis=1)
nvmax = np.max(nv[self.elements()],axis=1,keepdims=True)
h = self._get_matrix("element_size",self.n_elements(),1)
cfl = nvmax / (0.1*h)
......@@ -549,20 +566,24 @@ class FluidProblem :
f.restype = POINTER(typec)
return np.ctypeslib.as_array(f(self._fp),shape=(nrow,ncol))
def solution(self) :
"""Gives access to the fluid field value at the mesh nodes."""
sol = self.solution_flat()[:self.n_fields()*self.n_nodes()].reshape(-1, self.n_fields())
return sol
# def solution(self) :
# """Gives access to the fluid field value at the mesh nodes."""
# sol = self.solution_flat()[:self.n_fields()*self.n_nodes()].reshape(-1, self.n_fields())
# return sol
# return self._get_matrix("solution", self.n_nodes(), self.n_fields())
def solution_flat(self) :
def solution(self) :
"""Gives access to the fluid field value at each degree of freedom"""
ndof = np.max(self.set_numbering(self.order_fields))+1
return self._get_matrix("solution", ndof, 1).reshape([-1])
def velocity(self) :
"""Gives access to the fluid velocity value at the mesh nodes."""
return self.solution()[:,:-1]
"""Gives access to the fluid velocity value."""
ndof = np.max(self.set_numbering(self.order_fields))+1
id_v = np.r_[np.repeat(np.arange(self._dim)[None,:],self.n_nodes(),axis=0).reshape(-1)
+ self.n_fields()*np.repeat(np.arange(self.n_nodes()), self._dim),
np.arange(self.n_fields()*self.n_nodes(), ndof)]
return self.solution()[id_v].reshape(-1, self._dim)
def mesh_velocity(self) :
"""Give access to the mesh velocity value at the mesh nodes."""
......@@ -570,7 +591,8 @@ class FluidProblem :
def pressure(self) :
"""Gives access to the fluid field value at the mesh nodes."""
return self.solution()[:,[-1]]
id_p = np.arange(self._dim, self.n_fields()*self.n_nodes(), self.n_fields())
return self.solution()[id_p].reshape(-1)
def boundary_force(self) :
"""Give access to force exerted by the fluid on the boundaries."""
......@@ -679,57 +701,46 @@ class FluidProblem :
def local_map(self, idx):
self._lib.fluid_problem_local_map(self._fp, _np2c(idx, dtype=np.int32))
@timeit
# def set_numbering(self, order=[]):
# mesh_cl = MeshCl(self.elements())
# maps = [LagrangeSpace(mesh_cl, o).element_nodes for o in order]
# ndofs = [lmap.shape[1] for lmap in maps]
# ndofs_unique = np.unique(ndofs)
# n_order = len(order)
# # map en u0,v0,p0, u1,v1,p1, ...
# idx_uvp = np.full([self.n_elements(), np.sum(ndofs)], -1, dtype=np.int32)
# for i,ndof in enumerate(ndofs_unique):
# ndof_by_order = ndof if i == 0 else ndof-ndofs_unique[i-1]
# offset_dof = 0 if i == 0 else idx_uvp.max()+1
# offset_ind = np.where(idx_uvp[0]==-1)[0][0]
# order_to_use = np.where(ndofs >= ndof)[0]
# offset_order = 0 if i == 0 else ndof-ndofs_unique[i-1]
# n_order = len(order_to_use)
# for j in range(ndof_by_order):
# for k in range(n_order):
# lmap = maps[order_to_use[k]][:,offset_order:] - maps[order_to_use[k]][:,offset_order:].min()
# idx_uvp[:,offset_ind + j*n_order+k] = n_order*lmap[:,j]+k+offset_dof
# # map en u0,u1, ..., v0, v1, ..., p0, p1, ...
# idx = np.full([self.n_elements(), np.sum(ndofs)], -1, dtype=np.int32)
# n = ndofs_unique[0]
# for i in range(n_order):
# idx[:,i*n:(i+1)*n] = idx_uvp[:,i::n_order]
# return idx.reshape([-1])
def coordinates_p2(self):
connectivity = self.set_numbering([2]).reshape(self.n_elements(),-1)
x = np.zeros((connectivity.max()+1, self._dim))
x[connectivity[:,:self._dim+1]] = self.coordinates()[connectivity[:,:self._dim+1], :self._dim]
x[connectivity[:,self._dim+1]] = (x[connectivity[:,0]] + x[connectivity[:,1]])/2
x[connectivity[:,self._dim+2]] = (x[connectivity[:,1]] + x[connectivity[:,2]])/2
x[connectivity[:,self._dim+3]] = (x[connectivity[:,2]] + x[connectivity[:,0]])/2
return x if self._dim == 3 else np.insert(x, self._dim, 0, axis=1)
@timeit
def set_numbering(self, order=[]):
mesh_cl = MeshCl(self.elements())
maps = [LagrangeSpace(mesh_cl, o).element_nodes for o in order]
ndofs = [lmap.shape[1] for lmap in maps]
ndofs_unique = np.unique(ndofs)
n_order = len(order)
# map en u0,v0,p0, u1,v1,p1, ...
ndofs = np.array([lmap.shape[1] for lmap in maps])
dof_simplices = np.array([1, self._dim+1, 3*self._dim-3], dtype=np.int32)
nloc = np.array([
[1, 0, 0], #p0
[0, 1, 0], #p1
[0, 1, 1] #p2
], dtype=np.int32)
localsize = np.sum(nloc*dof_simplices[None,:], axis=1)
idx_uvp = np.full([self.n_elements(), np.sum(ndofs)], -1, dtype=np.int32)
for i,ndof in enumerate(ndofs_unique):
ndof_by_order = ndof if i == 0 else ndof-ndofs_unique[i-1]
offset_dof = 0 if i == 0 else idx_uvp.max()+1
offset_ind = np.where(idx_uvp[0]==-1)[0][0]
order_to_use = np.where(ndofs >= ndof)[0]
offset_order = 0 if i == 0 else ndof-ndofs_unique[i-1]
n_order = len(order_to_use)
for j in range(ndof_by_order):
for k in range(n_order):
lmap = maps[order_to_use[k]][:,offset_order:] - maps[order_to_use[k]][:,offset_order:].min()
idx_uvp[:,offset_ind + j*n_order+k] = n_order*lmap[:,j]+k+offset_dof
offsets = np.zeros(ndofs.shape[0], dtype=np.int32)
for s in range(0, dof_simplices.shape[0]):
has_simplices = np.where(nloc[order, s] != 0)[0]
if has_simplices.shape[0] == 0:
continue
n = has_simplices.shape[0]
map_dof = np.array([maps[j][:,offsets[has_simplices[j]]:offsets[has_simplices[j]]+localsize[s]] for j in has_simplices])
offsets = offsets[has_simplices]+localsize[s]
map_dof = np.moveaxis(map_dof, [0, 1, 2], [-1, 0, 1]).reshape(-1, n*dof_simplices[s]) - map_dof.min()
i = np.where(idx_uvp[0]==-1)[0].min()
end = i + dof_simplices[s]*n
idx_uvp[:,i:end] = np.repeat(np.arange(n)[None,:], dof_simplices[s], axis=0).reshape(-1)[None,:] \
+ n*map_dof + idx_uvp.max()+1
# map en u0,u1, ..., v0, v1, ..., p0, p1, ...
idx = np.full([self.n_elements(),np.sum(ndofs)],-1,dtype=np.int32)
if np.all(np.array(order)==1):
permutation = np.repeat(n_order*np.arange(ndofs[0])[None,:], n_order, axis=0).reshape([-1]) \
+ np.repeat(np.arange(n_order), ndofs[0])
if np.all(np.array(order)==order[0]):
permutation = np.repeat(ndofs.shape[0]*np.arange(ndofs[0])[None,:], ndofs.shape[0], axis=0).reshape([-1]) \
+ np.repeat(np.arange(ndofs.shape[0]), ndofs[0])
idx = idx_uvp[:,permutation]
if order==[2,2,1]: # 2D
permutation = np.array([0,3,6,9,11,13,1,4,7,10,12,14,2,5,8], dtype=np.int32)
......
......@@ -102,7 +102,7 @@ ii = 0
#
# FLUID PROBLEM
#
fluid = fluid.FluidProblem(2,g,[nu*rho],[rho], solver="petsc4py")
fluid = fluid.FluidProblem(2,g,[nu*rho],[rho], solver="petsc4py", solver_options="pc_type lu", order=[2,2,1])
# Set the mesh geometry for the fluid computation
fluid.load_msh("mesh.msh")
fluid.set_wall_boundary("Bottom")
......@@ -118,7 +118,11 @@ 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())
fluid.implicit_euler(dt, check_residual_norm=1e-3)
time_integration._advance_particles(p, fluid, dt, 5, 1e-8)
fluid.move_particles(p.position(), p.velocity(), p.contact_forces())
t += dt
# Output files writting
if ii %outf == 0 :
......
......@@ -68,14 +68,14 @@ class Poiseuille(unittest.TestCase) :
#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,nu*rho,rho,solver="petsc4py", order=[2,2,1])
fluid.set_mean_pressure(-100000)
fluid = mbfluid.FluidProblem(2,g,nu*rho,rho,solver="petsc4py", solver_options="pc_type lu", order=[2,2,1])
fluid.set_mean_pressure(10000)
fluid.load_msh("mesh.msh")
fluid.set_wall_boundary("Left",velocity=[0,0])
fluid.set_wall_boundary("Bottom",velocity=[0,0])
fluid.set_wall_boundary("Top",velocity=[vUp,0])
fluid.set_wall_boundary("Right",velocity=[0,0])
fluid.set_wall_boundary("Right",velocity=[0,0])
ii = 0
t = 0
......@@ -91,7 +91,7 @@ class Poiseuille(unittest.TestCase) :
tic = time.time()
while ii < 1000 :
#Fluid solver
time_integration.iterate(fluid,None,dt,check_residual_norm=1)
time_integration.iterate(fluid,None,dt,check_residual_norm=1e-1)
t += dt
#Output files writting
if ii %outf == 0 :
......@@ -100,8 +100,8 @@ class Poiseuille(unittest.TestCase) :
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
s = fluid.solution()
x = fluid.coordinates()
s = fluid.velocity()
x = fluid.coordinates_p2()
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]
......
......@@ -64,7 +64,7 @@ class Poiseuille(unittest.TestCase) :
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)
fluid = mbfluid.FluidProblem(2,g,nu*rho,rho, order=[2,2,1])
fluid = mbfluid.FluidProblem(2,g,nu*rho,rho, order=[2,2,1], solver_options="pc_type lu")
fluid.load_msh("mesh.msh")
fluid.set_open_boundary("LeftUp",velocity=[lambda x : 1/(20*mu)*x[:,1]*(1-x[:, 1]),0])
fluid.set_open_boundary("LeftDown",velocity=[lambda x : 1/(20*mu)*x[:,1]*(1-x[:, 1]),0])
......@@ -73,7 +73,6 @@ class Poiseuille(unittest.TestCase) :
fluid.set_open_boundary("RightUp",pressure=0)
fluid.set_open_boundary("RightDown",pressure=0)
ii = 0
t = 0
......@@ -82,12 +81,11 @@ class Poiseuille(unittest.TestCase) :
fluid.write_vtk(outputdir,0,0)
fluid.read_vtk(outputdir,0)
ii = 0
tic = time.perf_counter()
while ii < 100 :
#Fluid solver
time_integration.iterate(fluid,None,dt,check_residual_norm=1e-5)
time_integration.iterate(fluid,None,dt,check_residual_norm=1e-1)
t += dt
#Output files writting
if ii %outf == 0 :
......@@ -96,8 +94,8 @@ class Poiseuille(unittest.TestCase) :
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.perf_counter() - tic))
s = fluid.solution_reshape()
x = fluid.coordinates()
s = fluid.velocity()
x = fluid.coordinates_p2()
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)
......
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