xref: /libCEED/examples/python/ex3_volume.py (revision 02219a082eb38cf2d3edc97fbfe55fa395a4dc99) 
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 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