xref: /libCEED/python/ceed_operator.py (revision ed094490f53e580908aa80e9fe815a6fd76d7526) 
1# Copyright (c) 2017-2025, Lawrence Livermore National Security, LLC and other CEED contributors
2# All Rights Reserved. See the top-level LICENSE and NOTICE files for details.
3#
4# SPDX-License-Identifier: BSD-2-Clause
5#
6# This file is part of CEED:  http://github.com/ceed
7
8from _ceed_cffi import ffi, lib
9import tempfile
10from abc import ABC
11from .ceed_constants import REQUEST_IMMEDIATE, REQUEST_ORDERED, NOTRANSPOSE
12
13# ------------------------------------------------------------------------------
14
15
16class _OperatorBase(ABC):
17    """Ceed Operator: composed FE-type operations on vectors."""
18
19    # Destructor
20    def __del__(self):
21        # libCEED call
22        err_code = lib.CeedOperatorDestroy(self._pointer)
23        self._ceed._check_error(err_code)
24
25    # Representation
26    def __repr__(self):
27        return "<CeedOperator instance at " + hex(id(self)) + ">"
28
29    # String conversion for print() to stdout
30    def __str__(self):
31        """View an Operator via print()."""
32
33        # libCEED call
34        with tempfile.NamedTemporaryFile() as key_file:
35            with open(key_file.name, 'r+') as stream_file:
36                stream = ffi.cast("FILE *", stream_file)
37
38                err_code = lib.CeedOperatorView(self._pointer[0], stream)
39                self._ceed._check_error(err_code)
40
41                stream_file.seek(0)
42                out_string = stream_file.read()
43
44        return out_string
45
46    # Check Operator setup
47    def check(self):
48        """Check if a CeedOperator is ready to be used"""
49        # libCEED call
50        err_code = lib.CeedOperatorCheckReady(self._pointer[0])
51        self._ceed._check_error(err_code)
52
53    # Assemble linear diagonal
54    def linear_assemble_diagonal(self, d, request=REQUEST_IMMEDIATE):
55        """Assemble the diagonal of a square linear Operator
56
57           Args:
58             d: Vector to store assembled Operator diagonal
59             **request: Ceed request, default CEED_REQUEST_IMMEDIATE"""
60
61        # libCEED call
62        err_code = lib.CeedOperatorLinearAssembleDiagonal(self._pointer[0],
63                                                          d._pointer[0], request)
64        self._ceed._check_error(err_code)
65
66    # Assemble add linear diagonal
67    def linear_assemble_add_diagonal(self, d, request=REQUEST_IMMEDIATE):
68        """Sum the diagonal of a square linear Operator into a Vector
69
70           Args:
71             d: Vector to store assembled Operator diagonal
72             **request: Ceed request, default CEED_REQUEST_IMMEDIATE"""
73
74        # libCEED call
75        err_code = lib.CeedOperatorLinearAssembleAddDiagonal(self._pointer[0],
76                                                             d._pointer[0], request)
77        self._ceed._check_error(err_code)
78
79    # Assemble linear point block diagonal
80    def linear_assemble_point_block_diagonal(
81            self, d, request=REQUEST_IMMEDIATE):
82        """Assemble the point block diagonal of a square linear Operator
83
84           Args:
85             d: Vector to store assembled Operator point block diagonal,
86                  provided in row-major form with an ncomp*ncomp block
87                  at each node
88             **request: Ceed request, default CEED_REQUEST_IMMEDIATE"""
89
90        # libCEED call
91        err_code = lib.CeedOperatorLinearAssemblePointBlockDiagonal(self._pointer[0],
92                                                                    d._pointer[0], request)
93        self._ceed._check_error(err_code)
94
95    # Assemble linear point block diagonal
96    def linear_assemble_add_point_block_diagonal(
97            self, d, request=REQUEST_IMMEDIATE):
98        """Sum the point block diagonal of a square linear Operator into a Vector
99
100           Args:
101             d: Vector to store assembled Operator point block diagonal,
102                  provided in row-major form with an ncomp*ncomp block
103                  at each node
104             **request: Ceed request, default CEED_REQUEST_IMMEDIATE"""
105
106        # libCEED call
107        err_code = lib.CeedOperatorLinearAssembleAddPointBlockDiagonal(self._pointer[0],
108                                                                       d._pointer[0], request)
109        self._ceed._check_error(err_code)
110
111    # Set name
112    def name(self, name):
113        """Set name of Operator for print output
114
115           Args:
116             name: Name to set"""
117
118        name = ffi.new("char[]", name.encode('ascii'))
119        err_code = lib.CeedOperatorSetName(self._pointer[0], name)
120        self._ceed._check_error(err_code)
121
122    # Apply CeedOperator
123    def apply(self, u, v, request=REQUEST_IMMEDIATE):
124        """Apply Operator to a vector.
125
126           Args:
127             u: Vector containing input state or CEED_VECTOR_NONE if there are no
128                  active inputs
129             v: Vector to store result of applying operator (must be distinct from u)
130                  or CEED_VECTOR_NONE if there are no active outputs
131             **request: Ceed request, default CEED_REQUEST_IMMEDIATE"""
132
133        # libCEED call
134        err_code = lib.CeedOperatorApply(self._pointer[0], u._pointer[0], v._pointer[0],
135                                         request)
136        self._ceed._check_error(err_code)
137
138    # Apply CeedOperator
139    def apply_add(self, u, v, request=REQUEST_IMMEDIATE):
140        """Apply Operator to a vector and add result to output vector.
141
142           Args:
143             u: Vector containing input state or CEED_VECTOR_NONE if there are no
144                  active inputs
145             v: Vector to sum in result of applying operator (must be distinct from u)
146                  or CEED_VECTOR_NONE if there are no active outputs
147             **request: Ceed request, default CEED_REQUEST_IMMEDIATE"""
148
149        # libCEED call
150        err_code = lib.CeedOperatorApplyAdd(self._pointer[0], u._pointer[0], v._pointer[0],
151                                            request)
152        self._ceed._check_error(err_code)
153
154    # Create Multigrid Level
155    def multigrid_create(self, p_mult_fine, rstr_coarse, basis_coarse):
156        """ Create a multigrid coarse operator and level transfer operators
157           for a CeedOperator with a Lagrange tensor basis for the active basis
158
159           Args:
160             p_mult_fine: L-vector multiplicity in parallel gather/scatter
161             basis_coarse: Coarse grid active vector basis
162             degree_coarse: Coarse grid basis polynomial order"""
163
164        # Operator pointers
165        opCoarsePointer = ffi.new("CeedOperator *")
166        opProlongPointer = ffi.new("CeedOperator *")
167        opRestrictPointer = ffi.new("CeedOperator *")
168
169        # libCEED call
170        lib.CeedOperatorMultigridLevelCreate(self._pointer[0],
171                                             p_mult_fine._pointer[0],
172                                             rstr_coarse._pointer[0],
173                                             basis_coarse._pointer[0],
174                                             opCoarsePointer,
175                                             opProlongPointer,
176                                             opRestrictPointer)
177
178        # Wrap operators
179        opCoarse = _OperatorWrap(
180            self._ceed, opCoarsePointer)
181        opProlong = _OperatorWrap(
182            self._ceed, opProlongPointer)
183        opRestrict = _OperatorWrap(
184            self._ceed, opRestrictPointer)
185
186        # Return
187        return [opCoarse, opProlong, opRestrict]
188
189    # Create Multigrid Level
190    def multigrid_create_tensor_h1(self, p_mult_fine, rstr_coarse, basis_coarse,
191                                   interp_C_to_F):
192        """ Create a multigrid coarse operator and level transfer operators
193           for a CeedOperator with a non-tensor basis for the active basis
194
195           Args:
196             p_mult_fine: L-vector multiplicity in parallel gather/scatter
197             rstr_coarse: Coarse grid restriction
198             basis_coarse: Coarse grid active vector basis
199             interp_C_to_F: Matrix for coarse to fine interpolation"""
200
201       # Setup arguments
202        interpCtoF_pointer = ffi.new("CeedScalar *")
203        interpCtoF_pointer = ffi.cast(
204            "CeedScalar *",
205            interp_C_to_F.__array_interface__['data'][0])
206
207        # Operator pointers
208        opCoarsePointer = ffi.new("CeedOperator *")
209        opProlongPointer = ffi.new("CeedOperator *")
210        opRestrictPointer = ffi.new("CeedOperator *")
211
212        # libCEED call
213        lib.CeedOperatorMultigridLevelCreateTensorH1(self._pointer[0],
214                                                     p_mult_fine._pointer[0],
215                                                     rstr_coarse._pointer[0],
216                                                     basis_coarse._pointer[0],
217                                                     interpCtoF_pointer,
218                                                     opCoarsePointer,
219                                                     opProlongPointer,
220                                                     opRestrictPointer)
221
222        # Wrap operators
223        opCoarse = _OperatorWrap(
224            self._ceed, opCoarsePointer)
225        opProlong = _OperatorWrap(
226            self._ceed, opProlongPointer)
227        opRestrict = _OperatorWrap(
228            self._ceed, opRestrictPointer)
229
230        # Return
231        return [opCoarse, opProlong, opRestrict]
232
233    # Create Multigrid Level
234    def multigrid_create_h1(self, p_mult_fine, rstr_coarse, basis_coarse,
235                            interp_C_to_F):
236        """ Create a multigrid coarse operator and level transfer operators
237           for a CeedOperator with a Lagrange tensor basis for the active basis
238
239           Args:
240             p_mult_fine: L-vector multiplicity in parallel gather/scatter
241             rstr_coarse: Coarse grid restriction
242             basis_coarse: Coarse grid active vector basis
243             interp_C_to_F: Matrix for coarse to fine interpolation"""
244
245       # Setup arguments
246        interpCtoF_pointer = ffi.new("CeedScalar *")
247        interpCtoF_pointer = ffi.cast(
248            "CeedScalar *",
249            interp_C_to_F.__array_interface__['data'][0])
250
251        # Operator pointers
252        opCoarsePointer = ffi.new("CeedOperator *")
253        opProlongPointer = ffi.new("CeedOperator *")
254        opRestrictPointer = ffi.new("CeedOperator *")
255
256        # libCEED call
257        lib.CeedOperatorMultigridLevelCreateH1(self._pointer[0],
258                                               p_mult_fine._pointer[0],
259                                               rstr_coarse._pointer[0],
260                                               basis_coarse._pointer[0],
261                                               interpCtoF_pointer,
262                                               opCoarsePointer,
263                                               opProlongPointer,
264                                               opRestrictPointer)
265
266        # Wrap operators
267        opCoarse = _OperatorWrap(
268            self._ceed, opCoarsePointer)
269        opProlong = _OperatorWrap(
270            self._ceed, opProlongPointer)
271        opRestrict = _OperatorWrap(
272            self._ceed, opRestrictPointer)
273
274        # Return
275        return [opCoarse, opProlong, opRestrict]
276
277
278# ------------------------------------------------------------------------------
279
280
281class Operator(_OperatorBase):
282    """Ceed Operator: composed FE-type operations on vectors."""
283
284    # Constructor
285    def __init__(self, ceed, qf, dqf=None, dqfT=None):
286        # CeedOperator object
287        self._pointer = ffi.new("CeedOperator *")
288
289        # Reference to Ceed
290        self._ceed = ceed
291
292        # libCEED call
293        err_code = lib.CeedOperatorCreate(self._ceed._pointer[0], qf._pointer[0],
294                                          dqf._pointer[0] if dqf else ffi.NULL,
295                                          dqfT._pointer[0] if dqfT else ffi.NULL,
296                                          self._pointer)
297        self._ceed._check_error(err_code)
298
299    # Add field to CeedOperator
300    def set_field(self, fieldname, restriction, basis, vector):
301        """Provide a field to an Operator for use by its QFunction.
302
303           Args:
304             fieldname: name of the field (to be matched with the same name used
305                          by QFunction)
306             restriction: ElemRestriction
307             basis: Basis in which the field resides or CEED_BASIS_NONE
308                      if using CEED_EVAL_NONE
309             vector: Vector to be used by Operator or CEED_VECTOR_ACTIVE
310                       if field is active or CEED_VECTOR_NONE if using
311                       CEED_EVAL_WEIGHT in the QFunction"""
312
313        # libCEED call
314        fieldnameAscii = ffi.new("char[]", fieldname.encode('ascii'))
315        err_code = lib.CeedOperatorSetField(self._pointer[0], fieldnameAscii,
316                                            restriction._pointer[0], basis._pointer[0],
317                                            vector._pointer[0])
318        self._ceed._check_error(err_code)
319
320# ------------------------------------------------------------------------------
321
322
323class CompositeOperator(_OperatorBase):
324    """Ceed Composite Operator: composition of multiple Operators."""
325
326    # Constructor
327    def __init__(self, ceed):
328        # CeedOperator object
329        self._pointer = ffi.new("CeedOperator *")
330
331        # Reference to Ceed
332        self._ceed = ceed
333        # libCEED call
334        err_code = lib.CeedOperatorCreateComposite(
335            self._ceed._pointer[0], self._pointer)
336        self._ceed._check_error(err_code)
337
338    # Add sub operators
339    def add_sub(self, subop):
340        """Add a sub-operator to a composite CeedOperator.
341
342           Args:
343             subop: sub-operator Operator"""
344
345        # libCEED call
346        err_code = lib.CeedOperatorCompositeAddSub(
347            self._pointer[0], subop._pointer[0])
348        self._ceed._check_error(err_code)
349
350# ------------------------------------------------------------------------------
351
352
353class _OperatorWrap(Operator):
354    """Wrap a CeedOperator pointer in a Operator object."""
355
356    # Constructor
357    def __init__(self, ceed, pointer):
358        # CeedOperator object
359        self._pointer = pointer
360
361        # Reference to Ceed
362        self._ceed = ceed
363
364# ------------------------------------------------------------------------------
365