diff --git a/demo/lbb-stability/Makefile b/demo/lbb-stability/Makefile index 9a78350..f5f1d71 100644 --- a/demo/lbb-stability/Makefile +++ b/demo/lbb-stability/Makefile @@ -10,15 +10,14 @@ results-2d.json: lbb-stability.py results-3d.json: lbb-stability.py @for element in mini taylor-hood crouzeix-raviart; do \ for num_layers in 1 2 3; do \ - for vdegree in 1 2 3 4; do \ + for vdegree in 2 3 4; do \ python lbb-stability.py \ --dimension 3 \ --element $$element \ - --num-cells-min 4 \ - --num-cells-max 32 \ - --num-cells-step 2 \ + --log-num-cells-min 2 \ + --log-num-cells-max 8 \ --num-layers $$num_layers \ - --vdegree $$vdegree \ + --uvdegree $$vdegree \ --output $@ ; \ done \ done \ diff --git a/demo/lbb-stability/lbb-stability.py b/demo/lbb-stability/lbb-stability.py index 2942dbe..a146a2a 100644 --- a/demo/lbb-stability/lbb-stability.py +++ b/demo/lbb-stability/lbb-stability.py @@ -43,7 +43,7 @@ def compute_inf_sup_constant(spaces, b, inner_products): def run( - *, dimension, num_cells, element, num_layers=None, vdegree=None, vdegree_delta=None + *, dimension, num_cells, element, num_layers=None, uvdegree=None, pvdegree=None ): mesh = firedrake.UnitSquareMesh(num_cells, num_cells, diagonal="crossed") h = mesh.cell_sizes.dat.data_ro[:].min() @@ -62,8 +62,8 @@ def run( pressure_element = FiniteElement("DG", "triangle", 1) if dimension == 3: - cg_z = FiniteElement("CG", "interval", vdegree) - dg_z = FiniteElement("DG", "interval", vdegree - vdegree_delta) + cg_z = FiniteElement("CG", "interval", uvdegree) + dg_z = FiniteElement("DG", "interval", pvdegree) velocity_element = TensorProductElement(velocity_element, cg_z) pressure_element = TensorProductElement(pressure_element, dg_z) mesh = firedrake.ExtrudedMesh(mesh, num_layers) @@ -92,28 +92,25 @@ def run( parser.add_argument("--dimension", type=int, choices=[2, 3]) elements = ["taylor-hood", "crouzeix-raviart", "mini"] parser.add_argument("--element", choices=elements, default="taylor-hood") -parser.add_argument("--num-cells-min", type=int, default=4) -parser.add_argument("--num-cells-max", type=int, default=32) -parser.add_argument("--num-cells-step", type=int, default=1) +parser.add_argument("--log-num-cells-min", type=int, default=2) +parser.add_argument("--log-num-cells-max", type=int, default=6) parser.add_argument("--num-layers", type=int, default=1) -parser.add_argument("--vdegree", type=int, default=1) -parser.add_argument("--vdegree-delta", type=int, default=1) +parser.add_argument("--uvdegree", type=int, default=2) parser.add_argument("--output") args = parser.parse_args() stability_constants = [] -nmin, nmax, nstep = args.num_cells_min, args.num_cells_max, args.num_cells_step -print(f"Element: {args.element}") -print(f"#layers: {args.num_layers}") -print(f"vdegree: {args.vdegree}") -for num_cells in range(nmin, nmax + 1, nstep): +PETSc.Sys.Print(f"Element: {args.element}") +PETSc.Sys.Print(f"#layers: {args.num_layers}") +PETSc.Sys.Print(f"uvdegree: {args.uvdegree}") +for num_cells in 2 ** np.arange(args.log_num_cells_min, args.log_num_cells_max, 1): h, λ = run( dimension=args.dimension, num_cells=num_cells, element=args.element, num_layers=args.num_layers, - vdegree=args.vdegree, - vdegree_delta=args.vdegree_delta, + uvdegree=args.uvdegree, + pvdegree=args.uvdegree - 2, ) stability_constants.append((h, λ)) print(".", flush=True, end="") @@ -129,7 +126,13 @@ def run( "results": stability_constants, } if args.dimension == 3: - result.update({"num_layers": args.num_layers, "vdegree": args.vdegree}) + result.update( + { + "num_layers": args.num_layers, + "uvdegree": args.uvdegree, + "pvdegree": args.uvdegree - 2, + } + ) data.append(result) with open(args.output, "w") as output_file: json.dump(data, output_file)