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