1*5aed82e4SJeremy L Thompson // Copyright (c) 2017-2024, Lawrence Livermore National Security, LLC and other CEED contributors. 23d8e8822SJeremy L Thompson // All Rights Reserved. See the top-level LICENSE and NOTICE files for details. 3ed264d09SValeria Barra // 43d8e8822SJeremy L Thompson // SPDX-License-Identifier: BSD-2-Clause 5ed264d09SValeria Barra // 63d8e8822SJeremy L Thompson // This file is part of CEED: http://github.com/ceed 7ed264d09SValeria Barra 8ed264d09SValeria Barra /// @file 9ed264d09SValeria Barra /// libCEED QFunctions for mass operator example for a vector field on the sphere using PETSc 10ed264d09SValeria Barra 11f6b55d2cSvaleriabarra #ifndef bp2sphere_h 12f6b55d2cSvaleriabarra #define bp2sphere_h 13f6b55d2cSvaleriabarra 14c9c2c079SJeremy L Thompson #include <ceed.h> 15ed264d09SValeria Barra #include <math.h> 16ed264d09SValeria Barra 17e83e87a5Sjeremylt // ----------------------------------------------------------------------------- 18ed264d09SValeria Barra // This QFunction sets up the rhs and true solution for the problem 19ed264d09SValeria Barra // ----------------------------------------------------------------------------- 202b730f8bSJeremy L Thompson CEED_QFUNCTION(SetupMassRhs3)(void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 21ed264d09SValeria Barra // Inputs 229b072555Sjeremylt const CeedScalar *X = in[0], *q_data = in[1]; 23ed264d09SValeria Barra // Outputs 24ed264d09SValeria Barra CeedScalar *true_soln = out[0], *rhs = out[1]; 25ed264d09SValeria Barra 26ed264d09SValeria Barra // Context 27ed264d09SValeria Barra const CeedScalar *context = (const CeedScalar *)ctx; 28ed264d09SValeria Barra const CeedScalar R = context[0]; 29ed264d09SValeria Barra 30ed264d09SValeria Barra // Quadrature Point Loop 312b730f8bSJeremy L Thompson CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 32ed264d09SValeria Barra // Compute latitude 33ed264d09SValeria Barra const CeedScalar theta = asin(X[i + 2 * Q] / R); 34ed264d09SValeria Barra 359b072555Sjeremylt // Use absolute value of latitude for true solution 36ed264d09SValeria Barra // Component 1 37ed264d09SValeria Barra true_soln[i + 0 * Q] = fabs(theta); 38ed264d09SValeria Barra // Component 2 39ed264d09SValeria Barra true_soln[i + 1 * Q] = 2 * true_soln[i + 0 * Q]; 40ed264d09SValeria Barra // Component 3 41ed264d09SValeria Barra true_soln[i + 2 * Q] = 3 * true_soln[i + 0 * Q]; 42ed264d09SValeria Barra 43ed264d09SValeria Barra // Component 1 449b072555Sjeremylt rhs[i + 0 * Q] = q_data[i] * true_soln[i]; 45ed264d09SValeria Barra // Component 2 46ed264d09SValeria Barra rhs[i + 1 * Q] = 2 * rhs[i + 0 * Q]; 47ed264d09SValeria Barra // Component 3 48ed264d09SValeria Barra rhs[i + 2 * Q] = 3 * rhs[i + 0 * Q]; 49ed264d09SValeria Barra } // End of Quadrature Point Loop 50ed264d09SValeria Barra 51ed264d09SValeria Barra return 0; 52ed264d09SValeria Barra } 53ed264d09SValeria Barra 54e83e87a5Sjeremylt // ----------------------------------------------------------------------------- 55ed264d09SValeria Barra // This QFunction applies the mass operator for a vector field of 3 components. 56ed264d09SValeria Barra // 57ed264d09SValeria Barra // Inputs: 58ed264d09SValeria Barra // u - Input vector at quadrature points 599b072555Sjeremylt // q_data - Geometric factors 60ed264d09SValeria Barra // 61ed264d09SValeria Barra // Output: 62ed264d09SValeria Barra // v - Output vector (test functions) at quadrature points 63ed264d09SValeria Barra // ----------------------------------------------------------------------------- 642b730f8bSJeremy L Thompson CEED_QFUNCTION(Mass3)(void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 659b072555Sjeremylt const CeedScalar *u = in[0], *q_data = in[1]; 66ed264d09SValeria Barra CeedScalar *v = out[0]; 67ed264d09SValeria Barra 68ed264d09SValeria Barra // Quadrature Point Loop 692b730f8bSJeremy L Thompson CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 70ed264d09SValeria Barra // Component 1 719b072555Sjeremylt v[i + 0 * Q] = q_data[i] * u[i + 0 * Q]; 72ed264d09SValeria Barra // Component 2 739b072555Sjeremylt v[i + 1 * Q] = q_data[i] * u[i + 1 * Q]; 74ed264d09SValeria Barra // Component 3 759b072555Sjeremylt v[i + 2 * Q] = q_data[i] * u[i + 2 * Q]; 76ed264d09SValeria Barra } // End of Quadrature Point Loop 77ed264d09SValeria Barra 78ed264d09SValeria Barra return 0; 79ed264d09SValeria Barra } 80ed264d09SValeria Barra // ----------------------------------------------------------------------------- 81f6b55d2cSvaleriabarra 82f6b55d2cSvaleriabarra #endif // bp2sphere_h 83