|
| 1 | +#!/usr/bin/env python3 |
| 2 | +# Copyright (c) 2017-2025, Lawrence Livermore National Security, LLC and other CEED contributors. |
| 3 | +# All Rights Reserved. See the top-level LICENSE and NOTICE files for details. |
| 4 | +# |
| 5 | +# SPDX-License-Identifier: BSD-2-Clause |
| 6 | +# |
| 7 | +# This file is part of CEED: http://github.com/ceed |
| 8 | +# |
| 9 | +# libCEED example using diffusion operator to compute surface area |
| 10 | +# |
| 11 | +# Sample runs: |
| 12 | +# |
| 13 | +# python ex1_volume.py |
| 14 | +# python ex1_volume -c /cpu/self |
| 15 | +# python ex1_volume -c /gpu/cuda |
| 16 | + |
| 17 | +import sys |
| 18 | +import os |
| 19 | +import numpy as np |
| 20 | +import libceed |
| 21 | +import ex_common as common |
| 22 | + |
| 23 | + |
| 24 | +def main(): |
| 25 | + """Main function for volume example""" |
| 26 | + args = common.parse_arguments() |
| 27 | + return example_1(args) |
| 28 | + |
| 29 | + |
| 30 | +def example_1(args): |
| 31 | + """Compute volume using mass operator |
| 32 | +
|
| 33 | + Args: |
| 34 | + args: Parsed command line arguments |
| 35 | +
|
| 36 | + Returns: |
| 37 | + int: 0 on success, error code on failure |
| 38 | + """ |
| 39 | + # Process arguments |
| 40 | + dim = args.dim |
| 41 | + mesh_degree = max(args.mesh_degree, args.solution_degree) |
| 42 | + sol_degree = args.solution_degree |
| 43 | + num_qpts = args.quadrature_points |
| 44 | + problem_size = args.problem_size if args.problem_size > 0 else (8 * 16 if args.test else 256 * 1024) |
| 45 | + ncomp_x = dim # Number of coordinate components |
| 46 | + |
| 47 | + # Print configuration |
| 48 | + if not args.quiet: |
| 49 | + print("Selected options: [command line option] : <current value>") |
| 50 | + print(f" Ceed specification [-c] : {args.ceed}") |
| 51 | + print(f" Mesh dimension [-d] : {dim}") |
| 52 | + print(f" Mesh degree [-m] : {mesh_degree}") |
| 53 | + print(f" Solution degree [-p] : {sol_degree}") |
| 54 | + print(f" Num. 1D quadr. pts [-q] : {num_qpts}") |
| 55 | + print(f" Approx. # unknowns [-s] : {problem_size}") |
| 56 | + print(f" QFunction source [-g] : {'gallery' if args.gallery else 'user'}") |
| 57 | + |
| 58 | + # Initialize CEED |
| 59 | + ceed = libceed.Ceed(args.ceed) |
| 60 | + |
| 61 | + # Create bases |
| 62 | + # Tensor-product Lagrange basis for mesh coordinates |
| 63 | + mesh_basis = ceed.BasisTensorH1Lagrange( |
| 64 | + dim, ncomp_x, mesh_degree + 1, num_qpts, libceed.GAUSS) |
| 65 | + |
| 66 | + # Tensor-product Lagrange basis for solution |
| 67 | + solution_basis = ceed.BasisTensorH1Lagrange( |
| 68 | + dim, 1, sol_degree + 1, num_qpts, libceed.GAUSS) |
| 69 | + |
| 70 | + # Create mesh |
| 71 | + # Determine mesh size |
| 72 | + num_xyz = common.get_cartesian_mesh_size(dim, sol_degree, problem_size) |
| 73 | + if not args.quiet: |
| 74 | + print("\nMesh size : nx = %d" % num_xyz[0], end="") |
| 75 | + if dim > 1: |
| 76 | + print(", ny = %d" % num_xyz[1], end="") |
| 77 | + if dim > 2: |
| 78 | + print(", nz = %d" % num_xyz[2], end="") |
| 79 | + print() |
| 80 | + |
| 81 | + # Create element restrictions |
| 82 | + num_q_comp = 1 |
| 83 | + mesh_restriction, mesh_size, _, _, _ = common.build_cartesian_restriction( |
| 84 | + ceed, dim, num_xyz, mesh_degree, ncomp_x, num_q_comp, num_qpts, create_qdata=False) |
| 85 | + solution_restriction, sol_size, q_data_restriction, num_elem, elem_qpts = common.build_cartesian_restriction( |
| 86 | + ceed, dim, num_xyz, sol_degree, 1, num_q_comp, num_qpts, create_qdata=True) |
| 87 | + |
| 88 | + if not args.quiet: |
| 89 | + print("Number of mesh nodes : %d" % (mesh_size // dim)) |
| 90 | + print("Number of solution nodes : %d" % sol_size) |
| 91 | + |
| 92 | + # Create and transform mesh coordinates |
| 93 | + mesh_coords = ceed.Vector(mesh_size) |
| 94 | + common.set_cartesian_mesh_coords(ceed, dim, num_xyz, mesh_degree, mesh_coords) |
| 95 | + exact_volume, _ = common.transform_mesh_coords(dim, mesh_size, mesh_coords) |
| 96 | + |
| 97 | + # Create the QFunction that builds the mass operator (i.e. computes its quadrature data) and set its context data |
| 98 | + qf_build = None |
| 99 | + if args.gallery: |
| 100 | + qf_build = ceed.QFunctionByName(f"Mass{dim}DBuild") |
| 101 | + else: |
| 102 | + build_ctx = ceed.QFunctionContext() |
| 103 | + ctx_data = np.array([dim, dim], dtype=np.int32) |
| 104 | + build_ctx.set_data(ctx_data) |
| 105 | + |
| 106 | + qfs_so = common.load_qfs_so() |
| 107 | + file_dir = os.path.dirname(os.path.abspath(__file__)) |
| 108 | + |
| 109 | + qf_build = ceed.QFunction(1, qfs_so.build_mass, |
| 110 | + os.path.join(file_dir, "ex1-volume.h:build_mass")) |
| 111 | + qf_build.add_input("dx", dim * dim, libceed.EVAL_GRAD) |
| 112 | + qf_build.add_input("weights", 1, libceed.EVAL_WEIGHT) |
| 113 | + qf_build.add_output("qdata", num_q_comp, libceed.EVAL_NONE) |
| 114 | + qf_build.set_context(build_ctx) |
| 115 | + |
| 116 | + # Create the operator that builds the quadrature data for the mass operator |
| 117 | + op_build = ceed.Operator(qf_build) |
| 118 | + op_build.set_field("dx", mesh_restriction, mesh_basis, libceed.VECTOR_ACTIVE) |
| 119 | + op_build.set_field("weights", libceed.ELEMRESTRICTION_NONE, mesh_basis, libceed.VECTOR_NONE) |
| 120 | + op_build.set_field("qdata", q_data_restriction, libceed.BASIS_NONE, libceed.VECTOR_ACTIVE) |
| 121 | + |
| 122 | + # Compute the quadrature data for the mass operator |
| 123 | + q_data = ceed.Vector(num_elem * elem_qpts * num_q_comp) |
| 124 | + op_build.apply(mesh_coords, q_data) |
| 125 | + |
| 126 | + # Setup QFunction for applying the mass operator |
| 127 | + qf_mass = None |
| 128 | + if args.gallery: |
| 129 | + qf_mass = ceed.QFunctionByName("MassApply") |
| 130 | + else: |
| 131 | + build_ctx = ceed.QFunctionContext() |
| 132 | + ctx_data = np.array([dim, dim], dtype=np.int32) |
| 133 | + build_ctx.set_data(ctx_data) |
| 134 | + |
| 135 | + qfs_so = common.load_qfs_so() |
| 136 | + file_dir = os.path.dirname(os.path.abspath(__file__)) |
| 137 | + |
| 138 | + qf_mass = ceed.QFunction(1, qfs_so.apply_mass, |
| 139 | + os.path.join(file_dir, "ex1-volume.h:apply_mass")) |
| 140 | + qf_mass.add_input("u", 1, libceed.EVAL_INTERP) |
| 141 | + qf_mass.add_input("qdata", num_q_comp, libceed.EVAL_NONE) |
| 142 | + qf_mass.add_output("v", 1, libceed.EVAL_INTERP) |
| 143 | + qf_mass.set_context(build_ctx) |
| 144 | + |
| 145 | + # Create the mass operator |
| 146 | + op_mass = ceed.Operator(qf_mass) |
| 147 | + op_mass.set_field("u", solution_restriction, solution_basis, libceed.VECTOR_ACTIVE) |
| 148 | + op_mass.set_field("qdata", q_data_restriction, libceed.BASIS_NONE, q_data) |
| 149 | + op_mass.set_field("v", solution_restriction, solution_basis, libceed.VECTOR_ACTIVE) |
| 150 | + |
| 151 | + # Create solution vectors |
| 152 | + u = ceed.Vector(sol_size) |
| 153 | + v = ceed.Vector(sol_size) |
| 154 | + u.set_value(1.0) # Set all entries of u to 1.0 |
| 155 | + |
| 156 | + # Apply mass operator: v = M * u |
| 157 | + op_mass.apply(u, v) |
| 158 | + |
| 159 | + # Compute volume by summing all entries in v |
| 160 | + volume = 0.0 |
| 161 | + with v.array_read() as v_array: |
| 162 | + # Simply sum all values to compute the volume |
| 163 | + volume = np.sum(v_array) |
| 164 | + |
| 165 | + if not args.test: |
| 166 | + print() |
| 167 | + print(f"Exact mesh volume : {exact_volume:.14g}") |
| 168 | + print(f"Computed mesh volume : {volume:.14g}") |
| 169 | + print(f"Volume error : {volume - exact_volume:.14g}") |
| 170 | + else: |
| 171 | + # Test mode - check if error is within tolerance |
| 172 | + tol = 200 * libceed.EPSILON if dim == 1 else 1e-5 |
| 173 | + if abs(volume - exact_volume) > tol: |
| 174 | + print(f"Volume error : {volume - exact_volume:.14g}") |
| 175 | + sys.exit(1) |
| 176 | + |
| 177 | + return 0 |
| 178 | + |
| 179 | + |
| 180 | +if __name__ == "__main__": |
| 181 | + sys.exit(main()) |
0 commit comments