xref: /libCEED/tests/t505-operator.c (revision cae8b89ab4384ccb591c7143fd0ca8a192ac9d03) 
1*cae8b89aSjeremylt /// @file
2*cae8b89aSjeremylt /// Test CeedOperatorApplyAdd
3*cae8b89aSjeremylt /// \test Test CeedOperatorApplyAdd
4*cae8b89aSjeremylt #include <ceed.h>
5*cae8b89aSjeremylt #include <stdlib.h>
6*cae8b89aSjeremylt #include <math.h>
7*cae8b89aSjeremylt 
8*cae8b89aSjeremylt #include "t500-operator.h"
9*cae8b89aSjeremylt 
10*cae8b89aSjeremylt int main(int argc, char **argv) {
11*cae8b89aSjeremylt   Ceed ceed;
12*cae8b89aSjeremylt   CeedElemRestriction Erestrictx, Erestrictu, Erestrictxi, Erestrictui;
13*cae8b89aSjeremylt   CeedBasis bx, bu;
14*cae8b89aSjeremylt   CeedQFunction qf_setup, qf_mass;
15*cae8b89aSjeremylt   CeedOperator op_setup, op_mass;
16*cae8b89aSjeremylt   CeedVector qdata, X, U, V;
17*cae8b89aSjeremylt   const CeedScalar *hv;
18*cae8b89aSjeremylt   CeedInt nelem = 15, P = 5, Q = 8;
19*cae8b89aSjeremylt   CeedInt Nx = nelem+1, Nu = nelem*(P-1)+1;
20*cae8b89aSjeremylt   CeedInt indx[nelem*2], indu[nelem*P];
21*cae8b89aSjeremylt   CeedScalar x[Nx];
22*cae8b89aSjeremylt   CeedScalar sum;
23*cae8b89aSjeremylt 
24*cae8b89aSjeremylt   CeedInit(argv[1], &ceed);
25*cae8b89aSjeremylt 
26*cae8b89aSjeremylt   for (CeedInt i=0; i<Nx; i++)
27*cae8b89aSjeremylt     x[i] = (CeedScalar) i / (Nx - 1);
28*cae8b89aSjeremylt   for (CeedInt i=0; i<nelem; i++) {
29*cae8b89aSjeremylt     indx[2*i+0] = i;
30*cae8b89aSjeremylt     indx[2*i+1] = i+1;
31*cae8b89aSjeremylt   }
32*cae8b89aSjeremylt   // Restrictions
33*cae8b89aSjeremylt   CeedElemRestrictionCreate(ceed, nelem, 2, Nx, 1, CEED_MEM_HOST,
34*cae8b89aSjeremylt                             CEED_USE_POINTER, indx, &Erestrictx);
35*cae8b89aSjeremylt   CeedElemRestrictionCreateIdentity(ceed, nelem, 2, nelem*2, 1, &Erestrictxi);
36*cae8b89aSjeremylt 
37*cae8b89aSjeremylt   for (CeedInt i=0; i<nelem; i++) {
38*cae8b89aSjeremylt     for (CeedInt j=0; j<P; j++) {
39*cae8b89aSjeremylt       indu[P*i+j] = i*(P-1) + j;
40*cae8b89aSjeremylt     }
41*cae8b89aSjeremylt   }
42*cae8b89aSjeremylt   CeedElemRestrictionCreate(ceed, nelem, P, Nu, 1, CEED_MEM_HOST,
43*cae8b89aSjeremylt                             CEED_USE_POINTER, indu, &Erestrictu);
44*cae8b89aSjeremylt   CeedElemRestrictionCreateIdentity(ceed, nelem, Q, Q*nelem, 1, &Erestrictui);
45*cae8b89aSjeremylt 
46*cae8b89aSjeremylt   // Bases
47*cae8b89aSjeremylt   CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, 2, Q, CEED_GAUSS, &bx);
48*cae8b89aSjeremylt   CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, P, Q, CEED_GAUSS, &bu);
49*cae8b89aSjeremylt 
50*cae8b89aSjeremylt   // QFunctions
51*cae8b89aSjeremylt   CeedQFunctionCreateInterior(ceed, 1, setup, setup_loc, &qf_setup);
52*cae8b89aSjeremylt   CeedQFunctionAddInput(qf_setup, "_weight", 1, CEED_EVAL_WEIGHT);
53*cae8b89aSjeremylt   CeedQFunctionAddInput(qf_setup, "dx", 1, CEED_EVAL_GRAD);
54*cae8b89aSjeremylt   CeedQFunctionAddOutput(qf_setup, "rho", 1, CEED_EVAL_NONE);
55*cae8b89aSjeremylt 
56*cae8b89aSjeremylt   CeedQFunctionCreateInterior(ceed, 1, mass, mass_loc, &qf_mass);
57*cae8b89aSjeremylt   CeedQFunctionAddInput(qf_mass, "rho", 1, CEED_EVAL_NONE);
58*cae8b89aSjeremylt   CeedQFunctionAddInput(qf_mass, "u", 1, CEED_EVAL_INTERP);
59*cae8b89aSjeremylt   CeedQFunctionAddOutput(qf_mass, "v", 1, CEED_EVAL_INTERP);
60*cae8b89aSjeremylt 
61*cae8b89aSjeremylt   // Operators
62*cae8b89aSjeremylt   CeedOperatorCreate(ceed, qf_setup, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE,
63*cae8b89aSjeremylt                      &op_setup);
64*cae8b89aSjeremylt 
65*cae8b89aSjeremylt   CeedOperatorCreate(ceed, qf_mass, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE,
66*cae8b89aSjeremylt                      &op_mass);
67*cae8b89aSjeremylt 
68*cae8b89aSjeremylt   CeedVectorCreate(ceed, Nx, &X);
69*cae8b89aSjeremylt   CeedVectorSetArray(X, CEED_MEM_HOST, CEED_USE_POINTER, x);
70*cae8b89aSjeremylt   CeedVectorCreate(ceed, nelem*Q, &qdata);
71*cae8b89aSjeremylt 
72*cae8b89aSjeremylt   CeedOperatorSetField(op_setup, "_weight", Erestrictxi, CEED_NOTRANSPOSE,
73*cae8b89aSjeremylt                        bx, CEED_VECTOR_NONE);
74*cae8b89aSjeremylt   CeedOperatorSetField(op_setup, "dx", Erestrictx, CEED_NOTRANSPOSE,
75*cae8b89aSjeremylt                        bx, CEED_VECTOR_ACTIVE);
76*cae8b89aSjeremylt   CeedOperatorSetField(op_setup, "rho", Erestrictui, CEED_NOTRANSPOSE,
77*cae8b89aSjeremylt                        CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE);
78*cae8b89aSjeremylt 
79*cae8b89aSjeremylt   CeedOperatorSetField(op_mass, "rho", Erestrictui, CEED_NOTRANSPOSE,
80*cae8b89aSjeremylt                        CEED_BASIS_COLLOCATED, qdata);
81*cae8b89aSjeremylt   CeedOperatorSetField(op_mass, "u", Erestrictu, CEED_NOTRANSPOSE,
82*cae8b89aSjeremylt                        bu, CEED_VECTOR_ACTIVE);
83*cae8b89aSjeremylt   CeedOperatorSetField(op_mass, "v", Erestrictu, CEED_NOTRANSPOSE,
84*cae8b89aSjeremylt                        bu, CEED_VECTOR_ACTIVE);
85*cae8b89aSjeremylt 
86*cae8b89aSjeremylt   CeedOperatorApply(op_setup, X, qdata, CEED_REQUEST_IMMEDIATE);
87*cae8b89aSjeremylt 
88*cae8b89aSjeremylt   CeedVectorCreate(ceed, Nu, &U);
89*cae8b89aSjeremylt   CeedVectorSetValue(U, 1.0);
90*cae8b89aSjeremylt   CeedVectorCreate(ceed, Nu, &V);
91*cae8b89aSjeremylt   CeedVectorSetValue(V, 0.0);
92*cae8b89aSjeremylt 
93*cae8b89aSjeremylt   // Apply with V = 0
94*cae8b89aSjeremylt   CeedOperatorApplyAdd(op_mass, U, V, CEED_REQUEST_IMMEDIATE);
95*cae8b89aSjeremylt 
96*cae8b89aSjeremylt   // Check output
97*cae8b89aSjeremylt   CeedVectorGetArrayRead(V, CEED_MEM_HOST, &hv);
98*cae8b89aSjeremylt   sum = 0.;
99*cae8b89aSjeremylt   for (CeedInt i=0; i<Nu; i++)
100*cae8b89aSjeremylt     sum += hv[i];
101*cae8b89aSjeremylt   if (fabs(sum-1.)>1e-10) printf("Computed Area: %f != True Area: 1.0\n", sum);
102*cae8b89aSjeremylt   CeedVectorRestoreArrayRead(V, &hv);
103*cae8b89aSjeremylt 
104*cae8b89aSjeremylt   // Apply with V = 1
105*cae8b89aSjeremylt   CeedVectorSetValue(V, 1.0);
106*cae8b89aSjeremylt   CeedOperatorApplyAdd(op_mass, U, V, CEED_REQUEST_IMMEDIATE);
107*cae8b89aSjeremylt 
108*cae8b89aSjeremylt   // Check output
109*cae8b89aSjeremylt   CeedVectorGetArrayRead(V, CEED_MEM_HOST, &hv);
110*cae8b89aSjeremylt   sum = -Nu;
111*cae8b89aSjeremylt   for (CeedInt i=0; i<Nu; i++)
112*cae8b89aSjeremylt     sum += hv[i];
113*cae8b89aSjeremylt   if (fabs(sum-(1.))>1e-10) printf("Computed Area: %f != True Area: 1.0\n", sum);
114*cae8b89aSjeremylt   CeedVectorRestoreArrayRead(V, &hv);
115*cae8b89aSjeremylt 
116*cae8b89aSjeremylt   CeedQFunctionDestroy(&qf_setup);
117*cae8b89aSjeremylt   CeedQFunctionDestroy(&qf_mass);
118*cae8b89aSjeremylt   CeedOperatorDestroy(&op_setup);
119*cae8b89aSjeremylt   CeedOperatorDestroy(&op_mass);
120*cae8b89aSjeremylt   CeedElemRestrictionDestroy(&Erestrictu);
121*cae8b89aSjeremylt   CeedElemRestrictionDestroy(&Erestrictx);
122*cae8b89aSjeremylt   CeedElemRestrictionDestroy(&Erestrictui);
123*cae8b89aSjeremylt   CeedElemRestrictionDestroy(&Erestrictxi);
124*cae8b89aSjeremylt   CeedBasisDestroy(&bu);
125*cae8b89aSjeremylt   CeedBasisDestroy(&bx);
126*cae8b89aSjeremylt   CeedVectorDestroy(&X);
127*cae8b89aSjeremylt   CeedVectorDestroy(&U);
128*cae8b89aSjeremylt   CeedVectorDestroy(&V);
129*cae8b89aSjeremylt   CeedVectorDestroy(&qdata);
130*cae8b89aSjeremylt   CeedDestroy(&ceed);
131*cae8b89aSjeremylt   return 0;
132*cae8b89aSjeremylt }
133