xref: /libCEED/examples/python/ex1_volume.py (revision 5a526491291e2ef13670ec99232a2cb0069702e5) 
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
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    return example_1(args)
28
29
30def 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
180if __name__ == "__main__":
181    sys.exit(main())
182