1#!/usr/bin/env python3 2# Copyright (c) 2017-2026, 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 ex3_volume.py 14# python ex3_volume -c /cpu/self 15# python ex3_volume -c /gpu/cuda 16 17import sys 18import os 19import numpy as np 20import libceed 21import ex_common as common 22 23 24def main(): 25 """Main function for volume example""" 26 args = common.parse_arguments() 27 example_3(args) 28 29 30def example_3(args): 31 """Compute volume using mass and diff 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 # Check - Gallery not supported 59 if args.gallery: 60 print("Gallery QFunction not supported for example 3") 61 sys.exit(1) 62 63 # Initialize CEED 64 ceed = libceed.Ceed(args.ceed) 65 66 # Create bases 67 # Tensor-product Lagrange basis for mesh coordinates 68 mesh_basis = ceed.BasisTensorH1Lagrange( 69 dim, ncomp_x, mesh_degree + 1, num_qpts, libceed.GAUSS) 70 71 # Tensor-product Lagrange basis for solution 72 solution_basis = ceed.BasisTensorH1Lagrange( 73 dim, 1, sol_degree + 1, num_qpts, libceed.GAUSS) 74 75 # Create mesh 76 # Determine mesh size 77 num_xyz = common.get_cartesian_mesh_size(dim, sol_degree, problem_size) 78 if not args.quiet: 79 print("\nMesh size : nx = %d" % num_xyz[0], end="") 80 if dim > 1: 81 print(", ny = %d" % num_xyz[1], end="") 82 if dim > 2: 83 print(", nz = %d" % num_xyz[2], end="") 84 print() 85 86 # Create element restrictions 87 num_q_comp = 1 + dim * (dim + 1) // 2 88 mesh_restriction, mesh_size, _, _, _ = common.build_cartesian_restriction( 89 ceed, dim, num_xyz, mesh_degree, ncomp_x, num_q_comp, num_qpts, create_qdata=False) 90 solution_restriction, sol_size, q_data_restriction, num_elem, elem_qpts = common.build_cartesian_restriction( 91 ceed, dim, num_xyz, sol_degree, 1, num_q_comp, num_qpts, create_qdata=True) 92 93 if not args.quiet: 94 print("Number of mesh nodes : %d" % (mesh_size // dim)) 95 print("Number of solution nodes : %d" % sol_size) 96 97 # Create and transform mesh coordinates 98 mesh_coords = ceed.Vector(mesh_size) 99 common.set_cartesian_mesh_coords(ceed, dim, num_xyz, mesh_degree, mesh_coords) 100 exact_volume, _ = common.transform_mesh_coords(dim, mesh_size, mesh_coords) 101 102 # Create QFunction context 103 build_ctx = ceed.QFunctionContext() 104 ctx_data = np.array([dim, dim], dtype=np.int32) 105 build_ctx.set_data(ctx_data) 106 107 # Load QFunctions 108 qfs_so = common.load_qfs_so() 109 file_dir = os.path.dirname(os.path.abspath(__file__)) 110 111 # Create the QFunction that builds the mass + diffusion operator (i.e. 112 # computes its quadrature data) and set its context data 113 qf_build = ceed.QFunction(1, qfs_so.build_mass_diff, 114 os.path.join(file_dir, "ex3-volume.h:build_mass_diff")) 115 qf_build.add_input("dx", dim * dim, libceed.EVAL_GRAD) 116 qf_build.add_input("weights", 1, libceed.EVAL_WEIGHT) 117 qf_build.add_output("qdata", num_q_comp, libceed.EVAL_NONE) 118 qf_build.set_context(build_ctx) 119 120 # Create the operator that builds the quadrature data for the mass + diffusion operator 121 op_build = ceed.Operator(qf_build) 122 op_build.set_field("dx", mesh_restriction, mesh_basis, libceed.VECTOR_ACTIVE) 123 op_build.set_field("weights", libceed.ELEMRESTRICTION_NONE, mesh_basis, libceed.VECTOR_NONE) 124 op_build.set_field("qdata", q_data_restriction, libceed.BASIS_NONE, libceed.VECTOR_ACTIVE) 125 126 # Compute the quadrature data for the mass + diffusion operator 127 q_data = ceed.Vector(num_elem * elem_qpts * num_q_comp) 128 op_build.apply(mesh_coords, q_data) 129 130 # Create the QFunction that defines the action of the mass + diffusion operator 131 qf_apply = ceed.QFunction(1, qfs_so.apply_mass_diff, 132 os.path.join(file_dir, "ex3-volume.h:apply_mass_diff")) 133 qf_apply.add_input("u", 1, libceed.EVAL_INTERP) 134 qf_apply.add_input("du", dim, libceed.EVAL_GRAD) 135 qf_apply.add_input("qdata", num_q_comp, libceed.EVAL_NONE) 136 qf_apply.add_output("v", 1, libceed.EVAL_INTERP) 137 qf_apply.add_output("dv", dim, libceed.EVAL_GRAD) 138 qf_apply.set_context(build_ctx) 139 140 # Create the mass + diffusion operator 141 op_apply = ceed.Operator(qf_apply) 142 op_apply.set_field("u", solution_restriction, solution_basis, libceed.VECTOR_ACTIVE) 143 op_apply.set_field("du", solution_restriction, solution_basis, libceed.VECTOR_ACTIVE) 144 op_apply.set_field("qdata", q_data_restriction, libceed.BASIS_NONE, q_data) 145 op_apply.set_field("v", solution_restriction, solution_basis, libceed.VECTOR_ACTIVE) 146 op_apply.set_field("dv", solution_restriction, solution_basis, libceed.VECTOR_ACTIVE) 147 148 # Create solution vectors 149 u = ceed.Vector(sol_size) 150 v = ceed.Vector(sol_size) 151 u.set_value(1.0) # Set all entries of u to 1.0 152 153 # Apply mass + diffusion operator: v = (M + K) * u 154 op_apply.apply(u, v) 155 156 # Compute volume by summing all entries in v 157 volume = 0.0 158 with v.array_read() as v_array: 159 # Simply sum all values to compute the volume 160 volume = np.sum(v_array) 161 162 if not args.test: 163 print() 164 print(f"Exact mesh volume : {exact_volume:.14g}") 165 print(f"Computed mesh volume : {volume:.14g}") 166 print(f"Volume error : {volume - exact_volume:.14g}") 167 else: 168 # Test mode - check if error is within tolerance 169 tol = 200 * libceed.EPSILON if dim == 1 else 1e-5 170 if abs(volume - exact_volume) > tol: 171 print(f"Volume error : {volume - exact_volume:.14g}") 172 sys.exit(1) 173 174 return 0 175 176 177if __name__ == "__main__": 178 main() 179