xref: /libCEED/python/ceed_operator.py (revision a0154adecfab8547cdc0febbbf40ac009dbe9d1d) 
1# Copyright (c) 2017-2022, 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    # Apply CeedOperator
112    def apply(self, u, v, request=REQUEST_IMMEDIATE):
113        """Apply Operator to a vector.
114
115           Args:
116             u: Vector containing input state or CEED_VECTOR_NONE if there are no
117                  active inputs
118             v: Vector to store result of applying operator (must be distinct from u)
119                  or CEED_VECTOR_NONE if there are no active outputs
120             **request: Ceed request, default CEED_REQUEST_IMMEDIATE"""
121
122        # libCEED call
123        err_code = lib.CeedOperatorApply(self._pointer[0], u._pointer[0], v._pointer[0],
124                                         request)
125        self._ceed._check_error(err_code)
126
127    # Apply CeedOperator
128    def apply_add(self, u, v, request=REQUEST_IMMEDIATE):
129        """Apply Operator to a vector and add result to output vector.
130
131           Args:
132             u: Vector containing input state or CEED_VECTOR_NONE if there are no
133                  active inputs
134             v: Vector to sum in result of applying operator (must be distinct from u)
135                  or CEED_VECTOR_NONE if there are no active outputs
136             **request: Ceed request, default CEED_REQUEST_IMMEDIATE"""
137
138        # libCEED call
139        err_code = lib.CeedOperatorApplyAdd(self._pointer[0], u._pointer[0], v._pointer[0],
140                                            request)
141        self._ceed._check_error(err_code)
142
143    # Create Multigrid Level
144    def multigrid_create(self, p_mult_fine, rstr_coarse, basis_coarse):
145        """ Create a multigrid coarse operator and level transfer operators
146           for a CeedOperator with a Lagrange tensor basis for the active basis
147
148           Args:
149             p_mult_fine: L-vector multiplicity in parallel gather/scatter
150             basis_coarse: Coarse grid active vector basis
151             degree_coarse: Coarse grid basis polynomial order"""
152
153        # Operator pointers
154        opCoarsePointer = ffi.new("CeedOperator *")
155        opProlongPointer = ffi.new("CeedOperator *")
156        opRestrictPointer = ffi.new("CeedOperator *")
157
158        # libCEED call
159        lib.CeedOperatorMultigridLevelCreate(self._pointer[0],
160                                             p_mult_fine._pointer[0],
161                                             rstr_coarse._pointer[0],
162                                             basis_coarse._pointer[0],
163                                             opCoarsePointer,
164                                             opProlongPointer,
165                                             opRestrictPointer)
166
167        # Wrap operators
168        opCoarse = _OperatorWrap(
169            self._ceed, opCoarsePointer)
170        opProlong = _OperatorWrap(
171            self._ceed, opProlongPointer)
172        opRestrict = _OperatorWrap(
173            self._ceed, opRestrictPointer)
174
175        # Return
176        return [opCoarse, opProlong, opRestrict]
177
178    # Create Multigrid Level
179    def multigrid_create_tensor_h1(self, p_mult_fine, rstr_coarse, basis_coarse,
180                                   interp_C_to_F):
181        """ Create a multigrid coarse operator and level transfer operators
182           for a CeedOperator with a non-tensor basis for the active basis
183
184           Args:
185             p_mult_fine: L-vector multiplicity in parallel gather/scatter
186             rstr_coarse: Coarse grid restriction
187             basis_coarse: Coarse grid active vector basis
188             interp_C_to_F: Matrix for coarse to fine interpolation"""
189
190       # Setup arguments
191        interpCtoF_pointer = ffi.new("CeedScalar *")
192        interpCtoF_pointer = ffi.cast(
193            "CeedScalar *",
194            interp_C_to_F.__array_interface__['data'][0])
195
196        # Operator pointers
197        opCoarsePointer = ffi.new("CeedOperator *")
198        opProlongPointer = ffi.new("CeedOperator *")
199        opRestrictPointer = ffi.new("CeedOperator *")
200
201        # libCEED call
202        lib.CeedOperatorMultigridLevelCreateTensorH1(self._pointer[0],
203                                                     p_mult_fine._pointer[0],
204                                                     rstr_coarse._pointer[0],
205                                                     basis_coarse._pointer[0],
206                                                     interpCtoF_pointer,
207                                                     opCoarsePointer,
208                                                     opProlongPointer,
209                                                     opRestrictPointer)
210
211        # Wrap operators
212        opCoarse = _OperatorWrap(
213            self._ceed, opCoarsePointer)
214        opProlong = _OperatorWrap(
215            self._ceed, opProlongPointer)
216        opRestrict = _OperatorWrap(
217            self._ceed, opRestrictPointer)
218
219        # Return
220        return [opCoarse, opProlong, opRestrict]
221
222    # Create Multigrid Level
223    def multigrid_create_h1(self, p_mult_fine, rstr_coarse, basis_coarse,
224                            interp_C_to_F):
225        """ Create a multigrid coarse operator and level transfer operators
226           for a CeedOperator with a Lagrange tensor basis for the active basis
227
228           Args:
229             p_mult_fine: L-vector multiplicity in parallel gather/scatter
230             rstr_coarse: Coarse grid restriction
231             basis_coarse: Coarse grid active vector basis
232             interp_C_to_F: Matrix for coarse to fine interpolation"""
233
234       # Setup arguments
235        interpCtoF_pointer = ffi.new("CeedScalar *")
236        interpCtoF_pointer = ffi.cast(
237            "CeedScalar *",
238            interp_C_to_F.__array_interface__['data'][0])
239
240        # Operator pointers
241        opCoarsePointer = ffi.new("CeedOperator *")
242        opProlongPointer = ffi.new("CeedOperator *")
243        opRestrictPointer = ffi.new("CeedOperator *")
244
245        # libCEED call
246        lib.CeedOperatorMultigridLevelCreateH1(self._pointer[0],
247                                               p_mult_fine._pointer[0],
248                                               rstr_coarse._pointer[0],
249                                               basis_coarse._pointer[0],
250                                               interpCtoF_pointer,
251                                               opCoarsePointer,
252                                               opProlongPointer,
253                                               opRestrictPointer)
254
255        # Wrap operators
256        opCoarse = _OperatorWrap(
257            self._ceed, opCoarsePointer)
258        opProlong = _OperatorWrap(
259            self._ceed, opProlongPointer)
260        opRestrict = _OperatorWrap(
261            self._ceed, opRestrictPointer)
262
263        # Return
264        return [opCoarse, opProlong, opRestrict]
265
266
267# ------------------------------------------------------------------------------
268
269
270class Operator(_OperatorBase):
271    """Ceed Operator: composed FE-type operations on vectors."""
272
273    # Constructor
274    def __init__(self, ceed, qf, dqf=None, dqfT=None):
275        # CeedOperator object
276        self._pointer = ffi.new("CeedOperator *")
277
278        # Reference to Ceed
279        self._ceed = ceed
280
281        # libCEED call
282        err_code = lib.CeedOperatorCreate(self._ceed._pointer[0], qf._pointer[0],
283                                          dqf._pointer[0] if dqf else ffi.NULL,
284                                          dqfT._pointer[0] if dqfT else ffi.NULL,
285                                          self._pointer)
286        self._ceed._check_error(err_code)
287
288    # Add field to CeedOperator
289    def set_field(self, fieldname, restriction, basis, vector):
290        """Provide a field to an Operator for use by its QFunction.
291
292           Args:
293             fieldname: name of the field (to be matched with the same name used
294                          by QFunction)
295             restriction: ElemRestriction
296             basis: Basis in which the field resides or CEED_BASIS_COLLOCATED
297                      if collocated with quadrature points
298             vector: Vector to be used by Operator or CEED_VECTOR_ACTIVE
299                       if field is active or CEED_VECTOR_NONE if using
300                       CEED_EVAL_WEIGHT in the QFunction"""
301
302        # libCEED call
303        fieldnameAscii = ffi.new("char[]", fieldname.encode('ascii'))
304        err_code = lib.CeedOperatorSetField(self._pointer[0], fieldnameAscii,
305                                            restriction._pointer[0], basis._pointer[0],
306                                            vector._pointer[0])
307        self._ceed._check_error(err_code)
308
309# ------------------------------------------------------------------------------
310
311
312class CompositeOperator(_OperatorBase):
313    """Ceed Composite Operator: composition of multiple Operators."""
314
315    # Constructor
316    def __init__(self, ceed):
317        # CeedOperator object
318        self._pointer = ffi.new("CeedOperator *")
319
320        # Reference to Ceed
321        self._ceed = ceed
322        # libCEED call
323        err_code = lib.CeedCompositeOperatorCreate(
324            self._ceed._pointer[0], self._pointer)
325        self._ceed._check_error(err_code)
326
327    # Add sub operators
328    def add_sub(self, subop):
329        """Add a sub-operator to a composite CeedOperator.
330
331           Args:
332             subop: sub-operator Operator"""
333
334        # libCEED call
335        err_code = lib.CeedCompositeOperatorAddSub(
336            self._pointer[0], subop._pointer[0])
337        self._ceed._check_error(err_code)
338
339# ------------------------------------------------------------------------------
340
341
342class _OperatorWrap(Operator):
343    """Wrap a CeedOperator pointer in a Operator object."""
344
345    # Constructor
346    def __init__(self, ceed, pointer):
347        # CeedOperator object
348        self._pointer = pointer
349
350        # Reference to Ceed
351        self._ceed = ceed
352
353# ------------------------------------------------------------------------------
354