Skip to content

Commit

Permalink
Smaller mesh and cleanup
Browse files Browse the repository at this point in the history
  • Loading branch information
tlroy committed Dec 22, 2023
1 parent 44d2d90 commit 916af5f
Show file tree
Hide file tree
Showing 3 changed files with 104 additions and 85 deletions.
75 changes: 3 additions & 72 deletions examples/Navier_Stokes_irregular.py
Original file line number Diff line number Diff line change
@@ -1,23 +1,11 @@
# Navier-Stokes equations
# =======================
#

# TODO: These are nondimensionalized equations! We need the dimensional velocity

from firedrake import *
import sys
'''
N = 20
L=5
#M = RectangleMesh(N*L, N,L,0.2)
Ly = 0.1
Lx = 1.
#print(" Hauteur domaine=",pow(peclet,-1./3))
##mesh = RectangleBoundaryLayerMesh(50,50+int(10*pow(peclet,-1./3)),Lx,2*pow(peclet,-1./3)+Ly,50,1e-1, Ly_bdlayer = 5e-3, boundary=(3,1,))
##mesh = RectangleBoundaryLayerMesh(50,50,Lx,Ly,50,1e-1, Ly_bdlayer = 5e-3, boundary=(3,1,))
#M=RectangleMesh(50,50,Lx,Ly,quadrilateral=True)'''


from firedrake.petsc import PETSc


def navier_stokes_no_slip(mesh):
Expand Down Expand Up @@ -50,91 +38,34 @@ def navier_stokes_no_slip(mesh):
nullspace = MixedVectorSpaceBasis(
Z, [Z.sub(0), VectorSpaceBasis(constant=True)])

# Having set up the problem, we now move on to solving it. Some
# preconditioners, for example pressure convection-diffusion (PCD), require
# information about the the problem that is not easily accessible from
# the bilinear form. In the case of PCD, we need the Reynolds number
# and additionally which part of the mixed velocity-pressure space the
# velocity corresponds to. We provide this information to
# preconditioners by passing in a dictionary context to the solver.
# This is propagated down through the matrix-free operators and is
# therefore accessible to custom preconditioners. ::

appctx = {"Re": Re, "velocity_space": 0}

# Now we'll solve the problem. First, using a direct solver. Again, if
# MUMPS is not installed, this solve will not work, so we wrap the solve
# in a ``try/except`` block. ::

from firedrake.petsc import PETSc


# Now we'll show an example using the :class:`~.PCDPC` preconditioner
# that implements the pressure convection-diffusion approximation to the
# pressure Schur complement. We'll need more solver parameters this
# time, so again we'll set those up in a dictionary. ::

parameters = {"mat_type": "matfree",
"snes_monitor": None,

# We'll use a non-stationary Krylov solve for the Schur complement, so
# we need to use a flexible Krylov method on the outside. ::

"ksp_type": "fgmres",
"ksp_gmres_modifiedgramschmidt": None,
"ksp_monitor_true_residual": None,

# Now to configure the preconditioner::

"pc_type": "fieldsplit",
"pc_fieldsplit_type": "schur",
"pc_fieldsplit_schur_fact_type": "lower",

# we invert the velocity block with LU::

"fieldsplit_0_ksp_type": "preonly",
"fieldsplit_0_pc_type": "python",
"fieldsplit_0_pc_python_type": "firedrake.AssembledPC",
"fieldsplit_0_assembled_pc_type": "lu",

# and invert the schur complement inexactly using GMRES, preconditioned
# with PCD. ::

"fieldsplit_1_ksp_type": "gmres",
"fieldsplit_1_ksp_rtol": 1e-4,
"fieldsplit_1_pc_type": "python",
"fieldsplit_1_pc_python_type": "firedrake.PCDPC",

# We now need to configure the mass and stiffness solvers in the PCD
# preconditioner. For this example, we will just invert them with LU,
# although of course we can use a scalable method if we wish. First the
# mass solve::

"fieldsplit_1_pcd_Mp_ksp_type": "preonly",
"fieldsplit_1_pcd_Mp_pc_type": "lu",

# and the stiffness solve.::

"fieldsplit_1_pcd_Kp_ksp_type": "preonly",
"fieldsplit_1_pcd_Kp_pc_type": "lu",

# Finally, we just need to decide whether to apply the action of the
# pressure-space convection-diffusion operator with an assembled matrix
# or matrix free. Here we will use matrix-free::

"fieldsplit_1_pcd_Fp_mat_type": "matfree"}

# With the parameters set up, we can solve the problem, remembering to
# pass in the application context so that the PCD preconditioner can
# find the Reynolds number. ::

up.assign(0)

solve(F == 0, up, bcs=bcs, nullspace=nullspace, solver_parameters=parameters,
appctx=appctx)

# And finally we write the results to a file for visualisation. ::
#print(up.name)
u, p = up.subfunctions
u.rename("Velocity")
p.rename("Pressure")
Expand Down
17 changes: 4 additions & 13 deletions examples/Tworxn_irregular.py
Original file line number Diff line number Diff line change
Expand Up @@ -4,8 +4,6 @@
import sys
from Navier_Stokes_irregular import navier_stokes_no_slip


diffusion_coefficient=1.
peclet=10
damkohler=10
Ly = 0.1
Expand All @@ -23,7 +21,8 @@ def __init__(self):
Reactors. Industrial & Engineering Chemistry Research, 60(31),
pp.11824-11833.
"""
mesh=Mesh('squares.msh')
mesh=Mesh('squares_small.msh')
#mesh=Mesh('squares.msh')

C_1_inf = 1.
C_2_inf = Constant(0)
Expand Down Expand Up @@ -73,11 +72,6 @@ def set_velocity(self):
self.vel=navier_stokes_no_slip(self.mesh)





diffusion_coefficient=Lx/peclet
mass_transfert_coefficient=damkohler/Lx*diffusion_coefficient
solver = CarbonateSolver()
solver.setup_solver()
solver.solve()
Expand All @@ -87,14 +81,11 @@ def set_velocity(self):
flux = assemble(dot(grad(cC1), n) * ds(11))
flux1 = assemble(dot(grad(cC1), n)/dot(grad(cC1), n)*ds(11))

Peclet=Lx/diffusion_coefficient
Damkohler=mass_transfert_coefficient*Lx/diffusion_coefficient


print("Flux = %f" % flux)
print("surface : ",flux1)
fichier=open('results_square.dat','a')
fichier.write(str(Peclet)+' '+str(damkohler)+' '+str(flux)+' '+str(flux1)+' \n')
fichier.write(str(peclet)+' '+str(damkohler)+' '+str(flux)+' '+str(flux1)+' \n')
fichier.close()

File("SquareWave_"+str(Peclet)+str(damkohler)+".pvd").write(cC1)
File("SquareWave_"+str(peclet)+str(damkohler)+".pvd").write(cC1)
97 changes: 97 additions & 0 deletions examples/squares_small.geo
Original file line number Diff line number Diff line change
@@ -0,0 +1,97 @@
boundary_layer=5e-3;

size_boundary=0.1;
size_bas=boundary_layer/4;
size_bottom=boundary_layer/20;
height=0.3;
nb_square=2;
Amplitude=0.05;
//+
Point(1) = {1, 0, 0, size_bottom};
//+
Point(2) = {1, height, 0, size_boundary};
//+
Point(3) = {0, height, 0, size_boundary};
//+
Point(4) = {0, 0, 0, size_bottom};
Point(5) = {0,boundary_layer,0,size_bas};
Point(6) = {1,boundary_layer,0,size_bas};

k=5000;
For i In {1:nb_square}

Point(k) = {(1/4)/nb_square+(i-1)/nb_square+boundary_layer,boundary_layer,0,size_bas};
If(i>1)
Line(k-1)={k-1,k};
EndIf
Point(k+1) = {(1/4)/nb_square+(i-1)/nb_square+boundary_layer,boundary_layer-Amplitude,0,size_bas};
Point(k+2) = {(3/4)/nb_square+(i-1)/nb_square-boundary_layer,boundary_layer-Amplitude,0,size_bas};
Point(k+3) = {(3/4)/nb_square+(i-1)/nb_square-boundary_layer,boundary_layer,0,size_bas};
Line(k)={k,k+1};
Line(k+1)={k+1,k+2};
Line(k+2)={k+2,k+3};
k=k+4;
EndFor
Line(k-1)={k-1,6};
Line(4999)={5,5000};
max=k-1;


//+
Line(1) = {1, 6};
//+
Line(2) = {6, 2};
//+
Line(3) = {2, 3};
Line(4)={3,5};
Line(5)={5,4};
//+
k=7;


linep=7;
linec=1000;
list={5,6};

For i In {1:nb_square}
Point(k) = {(1/4)/nb_square+(i-1)/nb_square,0,0,size_bottom};
If(i>1)
Line(linec)={k-1,k};
list={list[],linec};
EndIf
Point(k+1) = {(1/4)/nb_square+(i-1)/nb_square,-Amplitude,0,size_bottom};
Point(k+2) = {(3/4)/nb_square+(i-1)/nb_square,-Amplitude,0,size_bottom};
Point(k+3) = {(3/4)/nb_square+(i-1)/nb_square,0,0,size_bottom};
Line(linec+1)={k,k+1};
Line(linep+1)={k+1,k+2};
Line(linec+2)={k+2,k+3};
list={list[],linec+1,linep+1,linec+2};
k=k+4;
linep=linep+1;
linec=linec+3;
EndFor
Line(6)={4,7};
Line(k-1)={k-1,1};
//+
Physical Curve("Cathode",11)={1001:linec-1};
Physical Curve("Pore",10)={8:linep};
Physical Curve("VerEdgeInlet", 12) = {4:5};
Physical Curve("VerEdgeOutlet", 14) = {1,2};
//+
Physical Curve("HorEdgeTop",13)={3};

Curve Loop(1)={2,3,4,4999:max};
//Curve Loop(2)={1,-max:-4999,5,6:k-1};
NumPoints = #list[];

Printf("The Array Contents are") ;
For index In {0:NumPoints-1}
Printf("%g ",list[index]) ;
EndFor
Curve Loop(2)={1,-max:-4999,list[],k-1};
Plane Surface(1)={1};
Plane Surface(2)={2};
Physical Surface(14)={1,2};

Mesh 2;

0 comments on commit 916af5f

Please sign in to comment.