diff --git a/examples/Hdiv-mixed/problems/richard2d.c b/examples/Hdiv-mixed/problems/richard2d.c index 89284b3d35..9033b0933b 100644 --- a/examples/Hdiv-mixed/problems/richard2d.c +++ b/examples/Hdiv-mixed/problems/richard2d.c @@ -19,16 +19,17 @@ #include "../include/register-problem.h" #include "../qfunctions/richard-system2d.h" +#include "../qfunctions/richard-force2d.h" #include "../qfunctions/pressure-boundary2d.h" PetscErrorCode Hdiv_RICHARD2D(Ceed ceed, ProblemData problem_data, void *ctx) { AppCtx app_ctx = *(AppCtx *)ctx; - //RICHARDContext richard_ctx; - //CeedQFunctionContext richard_context; + RICHARDContext richard_ctx; + CeedQFunctionContext richard_context; PetscFunctionBeginUser; - //PetscCall( PetscCalloc1(1, &richard_ctx) ); + PetscCall( PetscCalloc1(1, &richard_ctx) ); // ------------------------------------------------------ // SET UP POISSON_QUAD2D @@ -37,8 +38,8 @@ PetscErrorCode Hdiv_RICHARD2D(Ceed ceed, ProblemData problem_data, void *ctx) { problem_data->elem_node = 4; problem_data->q_data_size_face = 3; problem_data->quadrature_mode = CEED_GAUSS; - //problem_data->force = DarcyForce2D; - //problem_data->force_loc = DarcyForce2D_loc; + problem_data->force = RichardForce2D; + problem_data->force_loc = RichardForce2D_loc; problem_data->residual = RichardSystem2D; problem_data->residual_loc = RichardSystem2D_loc; problem_data->jacobian = JacobianRichardSystem2D; @@ -51,9 +52,38 @@ PetscErrorCode Hdiv_RICHARD2D(Ceed ceed, ProblemData problem_data, void *ctx) { // ------------------------------------------------------ // Command line Options // ------------------------------------------------------ + CeedScalar kappa = 10., alpha_a = 1., b_a = 10., rho_0 = 998.2, + beta = 0., g = 9.8, p0 = 101325; + PetscOptionsBegin(app_ctx->comm, NULL, "Options for Hdiv-mixed problem", NULL); + PetscCall( PetscOptionsScalar("-kappa", "Hydraulic Conductivity", NULL, + kappa, &kappa, NULL)); + PetscCall( PetscOptionsScalar("-alpha_a", "Parameter for relative permeability", + NULL, + alpha_a, &alpha_a, NULL)); + PetscCall( PetscOptionsScalar("-b_a", "Parameter for relative permeability", + NULL, + b_a, &b_a, NULL)); + PetscCall( PetscOptionsScalar("-rho_0", "Density at p0", NULL, + rho_0, &rho_0, NULL)); + PetscCall( PetscOptionsScalar("-beta", "Water compressibility", NULL, + beta, &beta, NULL)); PetscOptionsEnd(); + richard_ctx->kappa = kappa; + richard_ctx->alpha_a = alpha_a; + richard_ctx->b_a = b_a; + richard_ctx->rho_0 = rho_0; + richard_ctx->beta = beta; + richard_ctx->g = g; + richard_ctx->p0 = p0; + CeedQFunctionContextCreate(ceed, &richard_context); + CeedQFunctionContextSetData(richard_context, CEED_MEM_HOST, CEED_COPY_VALUES, + sizeof(*richard_ctx), richard_ctx); + problem_data->qfunction_context = richard_context; + CeedQFunctionContextSetDataDestroy(richard_context, CEED_MEM_HOST, + FreeContextPetsc); + //PetscCall( PetscFree(richard_ctx) ); PetscFunctionReturn(0); } diff --git a/examples/Hdiv-mixed/qfunctions/darcy-system2d.h b/examples/Hdiv-mixed/qfunctions/darcy-system2d.h index 4332c527cc..fc27f11d3f 100644 --- a/examples/Hdiv-mixed/qfunctions/darcy-system2d.h +++ b/examples/Hdiv-mixed/qfunctions/darcy-system2d.h @@ -184,4 +184,4 @@ CEED_QFUNCTION(JacobianDarcySystem2D)(void *ctx, CeedInt Q, // ----------------------------------------------------------------------------- -#endif //End of DARCY_MASS2D_H +#endif //End of DARCY_SYSTEM2D_H diff --git a/examples/Hdiv-mixed/qfunctions/darcy-system3d.h b/examples/Hdiv-mixed/qfunctions/darcy-system3d.h index 5def745c1f..4cd2adff55 100644 --- a/examples/Hdiv-mixed/qfunctions/darcy-system3d.h +++ b/examples/Hdiv-mixed/qfunctions/darcy-system3d.h @@ -195,4 +195,4 @@ CEED_QFUNCTION(JacobianDarcySystem3D)(void *ctx, CeedInt Q, // ----------------------------------------------------------------------------- -#endif //End of DARCY_MASS3D_H +#endif //End of DARCY_SYSTEM3D_H diff --git a/examples/Hdiv-mixed/qfunctions/richard-force2d.h b/examples/Hdiv-mixed/qfunctions/richard-force2d.h new file mode 100644 index 0000000000..32f651848e --- /dev/null +++ b/examples/Hdiv-mixed/qfunctions/richard-force2d.h @@ -0,0 +1,119 @@ +// Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at +// the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights +// reserved. See files LICENSE and NOTICE for details. +// +// This file is part of CEED, a collection of benchmarks, miniapps, software +// libraries and APIs for efficient high-order finite element and spectral +// element discretizations for exascale applications. For more information and +// source code availability see http://github.com/ceed. +// +// The CEED research is supported by the Exascale Computing Project 17-SC-20-SC, +// a collaborative effort of two U.S. Department of Energy organizations (Office +// of Science and the National Nuclear Security Administration) responsible for +// the planning and preparation of a capable exascale ecosystem, including +// software, applications, hardware, advanced system engineering and early +// testbed platforms, in support of the nation's exascale computing imperative. + +/// @file +/// Force of Richard problem 2D (quad element) using PETSc + +#ifndef RICHARD_FORCE2D_H +#define RICHARD_FORCE2D_H + +#include +#include "utils.h" + +// See Matthew Farthing, Christopher Kees, Cass Miller (2003) +// https://www.sciencedirect.com/science/article/pii/S0309170802001872 +// ----------------------------------------------------------------------------- +// Strong form: +// k*K^{-1} * u = -\grad(p) + rho*g in \Omega x [0,T] +// -\div(u) = -f + d (rho/rho_0*theta)/dt in \Omega x [0,T] +// p = p_b on \Gamma_D x [0,T] +// u.n = u_b on \Gamma_N x [0,T] +// p = p_0 in \Omega, t = 0 +// +// Where g is gravity vector, rho = rho_0*exp(beta * (p - p0)), p0 = 101325 Pa is atmospheric pressure +// f = fs/rho_0, where g is gravity, rho_0 is the density at p_0, K = kappa*I, and +// k_r = b_a + alpha_a * (\psi - x2), where \psi = p / (rho_0 * norm(g)) and x2 is vertical axis +// +// Weak form: Find (u, p) \in VxQ (V=H(div), Q=L^2) on \Omega +// (v, k*K^{-1} * u) -(v, rho*g) - (\div(v), p) = - _{\Gamma_D} +// -(q, \div(u)) = -(q, f) + (v, d (rho/rho_0*theta)/dt ) +// +// where k*K^{-1} = (rho_0^2*norm(g)/rho*k_r)*K^{-1} +// This QFunction sets up the force +// Inputs: +// x : interpolation of the physical coordinate +// w : weight of quadrature +// J : dx/dX. x physical coordinate, X reference coordinate [-1,1]^dim +// +// Output: +// force_u : which is 0.0 for this problem (- is in pressure-boundary qfunction) +// force_p : -(q, f) = -\int( q * f * w*detJ)dx +// ----------------------------------------------------------------------------- +// We have 3 experiment parameters as described in Table 1:P1, P2, P3 +// Matthew Farthing, Christopher Kees, Cass Miller (2003) +// https://www.sciencedirect.com/science/article/pii/S0309170802001872 +#ifndef RICHARD_CTX +#define RICHARD_CTX +typedef struct RICHARDContext_ *RICHARDContext; +struct RICHARDContext_ { + CeedScalar kappa; + CeedScalar alpha_a; + CeedScalar b_a; + CeedScalar rho_0; + CeedScalar beta; + CeedScalar g; + CeedScalar p0; +}; +#endif +// ----------------------------------------------------------------------------- +// Force evaluation for Richard problem +// ----------------------------------------------------------------------------- +CEED_QFUNCTION(RichardForce2D)(void *ctx, const CeedInt Q, + const CeedScalar *const *in, + CeedScalar *const *out) { + // *INDENT-OFF* + // Inputs + const CeedScalar (*coords) = in[0], + (*w) = in[1], + (*dxdX)[2][CEED_Q_VLA] = (const CeedScalar(*)[2][CEED_Q_VLA])in[2]; + // Outputs + CeedScalar (*rhs_u) = out[0], (*rhs_p) = out[1], + (*true_soln) = out[2]; + // Context + RICHARDContext context = (RICHARDContext)ctx; + const CeedScalar kappa = context->kappa;//10.; + // Quadrature Point Loop + CeedPragmaSIMD + for (CeedInt i=0; i +#include "utils.h" // See Matthew Farthing, Christopher Kees, Cass Miller (2003) // https://www.sciencedirect.com/science/article/pii/S0309170802001872 // ----------------------------------------------------------------------------- // Strong form: -// (rho_0^2*norm(g)/rho*k_r)*K^{-1} * u = -\grad(p) + rho*g in \Omega x [0,T] -// -\div(u) = -f + d (rho/rho_0*theta)/dt in \Omega x [0,T] -// p = p_b on \Gamma_D x [0,T] -// u.n = u_b on \Gamma_N x [0,T] -// p = p_0 in \Omega, t = 0 +// k*K^{-1} * u = -\grad(p) + rho*g in \Omega x [0,T] +// -\div(u) = -f + d (rho/rho_0*theta)/dt in \Omega x [0,T] +// p = p_b on \Gamma_D x [0,T] +// u.n = u_b on \Gamma_N x [0,T] +// p = p_0 in \Omega, t = 0 // -// Note: g is gravity vector, rho = rho_0*exp(beta * (p - p0)), p0 = 101325 Pa is atmospheric pressure -// f = fs/rho_0, where g is gravity, rho_0 is the density at p_0, K = K_star*I, and -// k_r = b_a + alpha_a * (psi - x2), where psi = p / (rho_0 * norm(g)) and x2 is vertical axis +// Where g is gravity vector, rho = rho_0*exp(beta * (p - p0)), p0 = 101325 Pa is atmospheric pressure +// f = fs/rho_0, where g is gravity, rho_0 is the density at p_0, K = kappa*I, and +// k_r = b_a + alpha_a * (\psi - x2), where \psi = p / (rho_0 * norm(g)) and x2 is vertical axis // // Weak form: Find (u, p) \in VxQ (V=H(div), Q=L^2) on \Omega -// (v, (rho_0^2*norm(g)/rho*k_r)*K^{-1} * u) - (\div(v), p) = (v, rho*g) - _{\Gamma_D} -// -(q, \div(u)) = -(q, f) + (v, d (rho*theta)/dt ) +// (v, k*K^{-1} * u) -(v, rho*g) - (\div(v), p) = - _{\Gamma_D} +// -(q, \div(u)) = -(q, f) + (v, d (rho/rho_0*theta)/dt ) // +// where k*K^{-1} = (rho_0^2*norm(g)/rho*k_r)*K^{-1} // This QFunction setup the mixed form of the above equation // Inputs: // w : weight of quadrature @@ -49,27 +51,31 @@ // p : basis_p at quadrature points // // Output: -// v : (v,u) = \int (v^T * u detJ*w) ==> \int (v^T J^T*J*u*w/detJ) +// v : (v, k*K^{-1} u) = \int (v^T * k*K^{-1}*u detJ*w) ==> \int (v^T J^T * k*K^{-1} *J*u*w/detJ) +// -(v, rho*g) = \int (v^T * rho*g detJ*w) ==> \int (v^T J^T * rho*g*w) // div(v) : -(\div(v), p) = -\int (div(v)^T * p *w) // q : -(q, \div(u)) = -\int (q^T * div(u) *w) // which create the following coupled system // D = [ M B^T ] // [ B 0 ] -// M = (v,u), B = -(q, \div(u)) -// Note we need to apply Piola map on the basis_u, which is J*u/detJ -// So (v,u) = \int (v^T * u detJ*w) ==> \int (v^T J^T*J*u*w/detJ) +// M = (v, k*K^{-1} u - rho*g), B = -(q, \div(u)) // ----------------------------------------------------------------------------- // We have 3 experiment parameters as described in Table 1:P1, P2, P3 // Matthew Farthing, Christopher Kees, Cass Miller (2003) // https://www.sciencedirect.com/science/article/pii/S0309170802001872 -typedef struct RICHARDP1Context_ *RICHARDP1Context; -struct RICHARDP1Context_ { - CeedScalar K_star; +#ifndef RICHARD_CTX +#define RICHARD_CTX +typedef struct RICHARDContext_ *RICHARDContext; +struct RICHARDContext_ { + CeedScalar kappa; CeedScalar alpha_a; CeedScalar b_a; CeedScalar rho_0; CeedScalar beta; + CeedScalar g; + CeedScalar p0; }; +#endif // ----------------------------------------------------------------------------- // Residual evaluation for Richard problem // ----------------------------------------------------------------------------- @@ -90,14 +96,14 @@ CEED_QFUNCTION(RichardSystem2D)(void *ctx, CeedInt Q, (*div_v) = (CeedScalar(*))out[1], (*q) = (CeedScalar(*))out[2]; // Context - //const RICHARDP1Context context = (RICHARDP1Context)ctx; - const CeedScalar K_star = 10.; - const CeedScalar alpha_a = 1.; - const CeedScalar b_a = 10.; - const CeedScalar rho_0 = 998.2; - const CeedScalar beta = 0.; - const CeedScalar g = 9.8; - const CeedScalar p0 = 101325; // atmospheric pressure + RICHARDContext context = (RICHARDContext)ctx; + const CeedScalar kappa = context->kappa;//10.; + const CeedScalar alpha_a = context->alpha_a;//1.; + const CeedScalar b_a = context->b_a;//10.; + const CeedScalar rho_0 = context->rho_0;//998.2; + const CeedScalar beta = context->beta;//0.; + const CeedScalar g = context->g;//9.8; + const CeedScalar p0 = context->p0;//101325; // atmospheric pressure // *INDENT-ON* // Quadrature Point Loop @@ -108,7 +114,7 @@ CEED_QFUNCTION(RichardSystem2D)(void *ctx, CeedInt Q, CeedScalar y = coords[i+1*Q]; const CeedScalar J[2][2] = {{dxdX[0][0][i], dxdX[1][0][i]}, {dxdX[0][1][i], dxdX[1][1][i]}}; - const CeedScalar detJ = J[0][0]*J[1][1] - J[0][1]*J[1][0]; + const CeedScalar det_J = MatDet2x2(J); // *INDENT-ON* // psi = p / (rho_0 * norm(g)) @@ -120,22 +126,32 @@ CEED_QFUNCTION(RichardSystem2D)(void *ctx, CeedInt Q, //k = rho_0^2*norm(g)/(rho*k_r) CeedScalar k = rho_0 * rho_0 * g / (rho * k_r); - // Piola map: J^T*k*K^{-1}*J*u*w/detJ - // Note K = K_star*I ==> K^{-1} = 1/K_star * I - // 1) Compute J^T*k*K^{-1}*J - CeedScalar JTkJ[2][2]; - for (CeedInt j = 0; j < 2; j++) { - for (CeedInt l = 0; l < 2; l++) { - JTkJ[j][l] = 0; - for (CeedInt m = 0; m < 2; m++) - JTkJ[j][l] += (k / K_star) * J[m][j] * J[m][l]; - } - } - // 2) Compute J^T*k*K^{-1}*J*u*w/detJ - for (CeedInt l = 0; l < 2; l++) { - v[l][i] = 0; - for (CeedInt m = 0; m < 2; m++) - v[l][i] += JTkJ[l][m] * u[m][i] * w[i]/detJ; + // (v, k*K^{-1} u): v = J^T*k*K^{-1}*J*u*w/detJ + // 1) Compute K^{-1}, note K = kappa*I + CeedScalar K[2][2] = {{kappa, 0.},{0., kappa}}; + const CeedScalar det_K = MatDet2x2(K); + CeedScalar K_inv[2][2]; + MatInverse2x2(K, det_K, K_inv); + + // 2) Compute k*K^{-1}*J + CeedScalar kKinv_J[2][2]; + AlphaMatMatMult2x2(k, K_inv, J, kKinv_J); + + // 3) Compute J^T * (k*K^{-1}*J) + CeedScalar JT_kKinv_J[2][2]; + AlphaMatTransposeMatMult2x2(1, J, kKinv_J, JT_kKinv_J); + + // 4) Compute (J^T*k*K^{-1}*J) * u * w /detJ + CeedScalar u1[2] = {u[0][i], u[1][i]}, v1[2]; + AlphaMatVecMult2x2(w[i]/det_J, JT_kKinv_J, u1, v1); + + // 5) -(v, rho*g): v = -J^T*rho*g*w + CeedScalar rho_g[2] = {0., rho*g}, v2[2]; + AlphaMatTransposeVecMult2x2(-w[i], J, rho_g, v2); + + // Output at quadrature points + for (CeedInt k = 0; k < 2; k++) { + v[k][i] = v1[k] + v2[k]; } div_v[i] = -p[i] * w[i]; @@ -167,14 +183,14 @@ CEED_QFUNCTION(JacobianRichardSystem2D)(void *ctx, CeedInt Q, (*div_dv) = (CeedScalar(*))out[1], (*dq) = (CeedScalar(*))out[2]; // Context - //const RICHARDP1Context context = (RICHARDP1Context)ctx; - const CeedScalar K_star = 10.; - const CeedScalar alpha_a = 1.; - const CeedScalar b_a = 10.; - const CeedScalar rho_0 = 998.2; - const CeedScalar beta = 0.; - const CeedScalar g = 9.8; - const CeedScalar p0 = 101325; // atmospheric pressure + RICHARDContext context = (RICHARDContext)ctx; + const CeedScalar kappa = context->kappa; + const CeedScalar alpha_a = context->alpha_a; + const CeedScalar b_a = context->b_a; + const CeedScalar rho_0 = context->rho_0; + const CeedScalar beta = context->beta; + const CeedScalar g = context->g; + const CeedScalar p0 = context->p0;// atmospheric pressure // *INDENT-ON* // Quadrature Point Loop @@ -185,7 +201,7 @@ CEED_QFUNCTION(JacobianRichardSystem2D)(void *ctx, CeedInt Q, CeedScalar y = coords[i+1*Q]; const CeedScalar J[2][2] = {{dxdX[0][0][i], dxdX[1][0][i]}, {dxdX[0][1][i], dxdX[1][1][i]}}; - const CeedScalar detJ = J[0][0]*J[1][1] - J[0][1]*J[1][0]; + const CeedScalar det_J = MatDet2x2(J); // *INDENT-ON* // psi = p / (rho_0 * norm(g)) @@ -198,27 +214,42 @@ CEED_QFUNCTION(JacobianRichardSystem2D)(void *ctx, CeedInt Q, CeedScalar k = rho_0 * rho_0 * g / (rho * k_r); // Piola map: J^T*k*K^{-1}*J*u*w/detJ - // Note K = K_star*I ==> K^{-1} = 1/K_star * I // The jacobian term - // dv = J^T* (k*K^{-1}) *J*du*w/detJ - J^T*(k*K^{-1} [(rho*k_r)/dp]*dp/(rho*k_r)) *J*u - - // 1) Compute J^T* (k*K^{-1}) *J - CeedScalar JTkJ[2][2]; - for (CeedInt j = 0; j < 2; j++) { - for (CeedInt l = 0; l < 2; l++) { - JTkJ[j][l] = 0; - for (CeedInt m = 0; m < 2; m++) - JTkJ[j][l] += (k / K_star) * J[m][j] * J[m][l]; - } - } + // dv = J^T* (k*K^{-1}) *J*du*w/detJ - [(rho*k_r),p*dp/(rho*k_r)]*J^T*(k*K^{-1}) *J*u*w/detJ + // -J^T * (beta*rho*g)*dp + // 1) Compute K^{-1}, note K = kappa*I + CeedScalar K[2][2] = {{kappa, 0.},{0., kappa}}; + const CeedScalar det_K = MatDet2x2(K); + CeedScalar K_inv[2][2]; + MatInverse2x2(K, det_K, K_inv); + + // 2) Compute k*K^{-1}*J + CeedScalar kKinv_J[2][2]; + AlphaMatMatMult2x2(k, K_inv, J, kKinv_J); + + // 3) Compute J^T * (k*K^{-1}*J) + CeedScalar JT_kKinv_J[2][2]; + AlphaMatTransposeMatMult2x2(1, J, kKinv_J, JT_kKinv_J); + + // 4) Compute (J^T*k*K^{-1}*J) * du * w /detJ + CeedScalar du1[2] = {du[0][i], du[1][i]}, dv1[2]; + AlphaMatVecMult2x2(w[i]/det_J, JT_kKinv_J, du1, dv1); + + // 5) Compute -(rho*k_r),p*dp/(rho*k_r)) + // (rho*k_r),p*dp = beta*rho*dp*k_r + rho*alpha*dp/(rho_0*norm(g)) + CeedScalar d_rhokr_dp = -(beta + alpha_a/(rho_0*g*k_r))*dp[i]; + + // 6) -[(rho*k_r),p*dp/(rho*k_r)]*J^T*(k*K^{-1}) *J*u*w/detJ + CeedScalar u1[2] = {u[0][i], u[1][i]}, dv2[2]; + AlphaMatVecMult2x2((d_rhokr_dp*w[i])/det_J, JT_kKinv_J, u1, dv2); + + // 7) -(v, rho*g): dv = -J^T * (beta*rho*g*dp)*w + CeedScalar drho_g_dp[2] = {0., beta *rho *g *dp[i]}, dv3[2]; + AlphaMatTransposeVecMult2x2(-w[i], J, drho_g_dp, dv3); - // 2) Compute [(rho*k_r)/dp]*dp/(rho*k_r)) - CeedScalar d_rhokr_dp = (beta + alpha_a/(rho_0*g*k_r))*dp[i]; - // 3) Compute dv - for (CeedInt l = 0; l < 2; l++) { - dv[l][i] = 0; - for (CeedInt m = 0; m < 2; m++) - dv[l][i] += JTkJ[l][m] * du[m][i] * w[i]/detJ - JTkJ[l][m]*d_rhokr_dp*u[m][i]; + // Output at quadrature points + for (CeedInt k = 0; k < 2; k++) { + dv[k][i] = dv1[k] + dv2[k] + dv3[k]; } div_dv[i] = -dp[i] * w[i]; @@ -230,4 +261,4 @@ CEED_QFUNCTION(JacobianRichardSystem2D)(void *ctx, CeedInt Q, // ----------------------------------------------------------------------------- -#endif //End of DARCY_MASS2D_H +#endif //End of RICHARD_SYSTEM2D_H diff --git a/examples/Hdiv-mixed/qfunctions/utils.h b/examples/Hdiv-mixed/qfunctions/utils.h index 449b509720..d19662fbc4 100644 --- a/examples/Hdiv-mixed/qfunctions/utils.h +++ b/examples/Hdiv-mixed/qfunctions/utils.h @@ -172,4 +172,19 @@ CEED_QFUNCTION_HELPER int AlphaMatVecMult2x2(const CeedScalar alpha, return 0; }; +// ----------------------------------------------------------------------------- +// Compute matrix-vector product: alpha*A^T*u +// ----------------------------------------------------------------------------- +CEED_QFUNCTION_HELPER int AlphaMatTransposeVecMult2x2(const CeedScalar alpha, + const CeedScalar A[2][2], const CeedScalar u[2], CeedScalar v[2]) { + // Compute v = alpha*A*u + for (CeedInt k = 0; k < 2; k++) { + v[k] = 0; + for (CeedInt m = 0; m < 2; m++) + v[k] += A[m][k] * u[m] * alpha; + } + + return 0; +}; + #endif // utils_qf_h diff --git a/interface/ceed-operator.c b/interface/ceed-operator.c index 3b8de1cead..02b730fe92 100644 --- a/interface/ceed-operator.c +++ b/interface/ceed-operator.c @@ -88,7 +88,6 @@ static int CeedOperatorCheckField(Ceed ceed, CeedQFunctionField qf_field, // LCOV_EXCL_STOP } - CeedFESpace FE_space = b->basis_space; // Field size switch(eval_mode) { case CEED_EVAL_NONE: