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