1*9ba83ac0SJeremy L Thompson // Copyright (c) 2017-2026, 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 11c0b5abf0SJeremy L Thompson #include <ceed/types.h> 12c0b5abf0SJeremy L Thompson #ifndef CEED_RUNNING_JIT_PASS 13ed264d09SValeria Barra #include <math.h> 14c0b5abf0SJeremy L Thompson #endif 15ed264d09SValeria Barra 16e83e87a5Sjeremylt // ----------------------------------------------------------------------------- 17ed264d09SValeria Barra // This QFunction sets up the rhs and true solution for the problem 18ed264d09SValeria Barra // ----------------------------------------------------------------------------- 192b730f8bSJeremy L Thompson CEED_QFUNCTION(SetupMassRhs3)(void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 20ed264d09SValeria Barra // Inputs 219b072555Sjeremylt const CeedScalar *X = in[0], *q_data = in[1]; 22ed264d09SValeria Barra // Outputs 23ed264d09SValeria Barra CeedScalar *true_soln = out[0], *rhs = out[1]; 24ed264d09SValeria Barra 25ed264d09SValeria Barra // Context 26ed264d09SValeria Barra const CeedScalar *context = (const CeedScalar *)ctx; 27ed264d09SValeria Barra const CeedScalar R = context[0]; 28ed264d09SValeria Barra 29ed264d09SValeria Barra // Quadrature Point Loop 302b730f8bSJeremy L Thompson CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 31ed264d09SValeria Barra // Compute latitude 32ed264d09SValeria Barra const CeedScalar theta = asin(X[i + 2 * Q] / R); 33ed264d09SValeria Barra 349b072555Sjeremylt // Use absolute value of latitude for true solution 35ed264d09SValeria Barra // Component 1 36ed264d09SValeria Barra true_soln[i + 0 * Q] = fabs(theta); 37ed264d09SValeria Barra // Component 2 38ed264d09SValeria Barra true_soln[i + 1 * Q] = 2 * true_soln[i + 0 * Q]; 39ed264d09SValeria Barra // Component 3 40ed264d09SValeria Barra true_soln[i + 2 * Q] = 3 * true_soln[i + 0 * Q]; 41ed264d09SValeria Barra 42ed264d09SValeria Barra // Component 1 439b072555Sjeremylt rhs[i + 0 * Q] = q_data[i] * true_soln[i]; 44ed264d09SValeria Barra // Component 2 45ed264d09SValeria Barra rhs[i + 1 * Q] = 2 * rhs[i + 0 * Q]; 46ed264d09SValeria Barra // Component 3 47ed264d09SValeria Barra rhs[i + 2 * Q] = 3 * rhs[i + 0 * Q]; 48ed264d09SValeria Barra } // End of Quadrature Point Loop 49ed264d09SValeria Barra 50ed264d09SValeria Barra return 0; 51ed264d09SValeria Barra } 52ed264d09SValeria Barra 53e83e87a5Sjeremylt // ----------------------------------------------------------------------------- 54ed264d09SValeria Barra // This QFunction applies the mass operator for a vector field of 3 components. 55ed264d09SValeria Barra // 56ed264d09SValeria Barra // Inputs: 57ed264d09SValeria Barra // u - Input vector at quadrature points 589b072555Sjeremylt // q_data - Geometric factors 59ed264d09SValeria Barra // 60ed264d09SValeria Barra // Output: 61ed264d09SValeria Barra // v - Output vector (test functions) at quadrature points 62ed264d09SValeria Barra // ----------------------------------------------------------------------------- 632b730f8bSJeremy L Thompson CEED_QFUNCTION(Mass3)(void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 649b072555Sjeremylt const CeedScalar *u = in[0], *q_data = in[1]; 65ed264d09SValeria Barra CeedScalar *v = out[0]; 66ed264d09SValeria Barra 67ed264d09SValeria Barra // Quadrature Point Loop 682b730f8bSJeremy L Thompson CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 69ed264d09SValeria Barra // Component 1 709b072555Sjeremylt v[i + 0 * Q] = q_data[i] * u[i + 0 * Q]; 71ed264d09SValeria Barra // Component 2 729b072555Sjeremylt v[i + 1 * Q] = q_data[i] * u[i + 1 * Q]; 73ed264d09SValeria Barra // Component 3 749b072555Sjeremylt v[i + 2 * Q] = q_data[i] * u[i + 2 * Q]; 75ed264d09SValeria Barra } // End of Quadrature Point Loop 76ed264d09SValeria Barra 77ed264d09SValeria Barra return 0; 78ed264d09SValeria Barra } 79ed264d09SValeria Barra // ----------------------------------------------------------------------------- 80