1!----------------------------------------------------------------------- 2 subroutine setup_mass(ctx,q,u1,u2,u3,u4,u5,u6,u7,u8,u9,u10,u11,u12,u13,& 3& u14,u15,u16,v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,v13,v14,v15,& 4& v16,ierr) 5 real*8 ctx 6 real*8 u1(1) 7 real*8 u2(1) 8 real*8 v1(1) 9 integer q,ierr 10 11 do i=1,q 12 v1(i)=u2(i)*(u1(i+q*0)*u1(i+q*3)-u1(i+q*1)*u1(i+q*2)) 13 enddo 14 15 ierr=0 16 end 17!----------------------------------------------------------------------- 18 subroutine setup_diff(ctx,q,u1,u2,u3,u4,u5,u6,u7,u8,u9,u10,u11,u12,u13,& 19& u14,u15,u16,v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,v13,v14,v15,& 20& v16,ierr) 21 real*8 ctx 22 real*8 u1(1) 23 real*8 u2(1) 24 real*8 v1(1) 25 real*8 w 26 integer q,ierr 27 28 do i=1,q 29 w=u2(i)/(u1(i+q*0)*u1(i+q*3)-u1(i+q*1)*u1(i+q*2)) 30 v1(i+q*0)=w*(u1(i+q*2)*u1(i+q*2)+u1(i+q*3)*u1(i+q*3)) 31 v1(i+q*1)=w*(u1(i+q*0)*u1(i+q*0)+u1(i+q*1)*u1(i+q*1)) 32 v1(i+q*2)=-w*(u1(i+q*0)*u1(i+q*2)+u1(i+q*2)*u1(i+q*3)) 33 enddo 34 35 ierr=0 36 end 37!----------------------------------------------------------------------- 38 subroutine apply(ctx,q,u1,u2,u3,u4,u5,u6,u7,u8,u9,u10,u11,u12,u13,u14,& 39& u15,u16,v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,v13,v14,v15,v16,ierr) 40 real*8 ctx 41 real*8 u1(1) 42 real*8 u2(1) 43 real*8 u3(1) 44 real*8 u4(1) 45 real*8 v1(1) 46 real*8 v2(1) 47 real*8 du0,du1 48 integer q,ierr 49 50 do i=1,q 51! mass 52 v1(i) = u2(i)*u4(i) 53! diff 54 du0=u1(i+q*0) 55 du1=u1(i+q*1) 56 v2(i+q*0)=u3(i+q*0)*du0+u3(i+q*2)*du1 57 v2(i+q*1)=u3(i+q*2)*du0+u3(i+q*1)*du1 58 enddo 59 60 ierr=0 61 end 62!----------------------------------------------------------------------- 63 subroutine apply_lin(ctx,q,u1,u2,u3,u4,u5,u6,u7,u8,u9,u10,u11,u12,u13,& 64& u14,u15,u16,v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,v13,v14,v15,& 65& v16,ierr) 66 real*8 ctx 67 real*8 u1(1) 68 real*8 u2(1) 69 real*8 u3(1) 70 real*8 v1(1) 71 real*8 v2(1) 72 real*8 du0,du1 73 integer q,ierr 74 75 do i=1,q 76 du0=u1(i+q*0) 77 du1=u1(i+q*1) 78 v1(i+q*0)=u2(i+q*0)*du0+u2(i+q*3)*du1+u2(i+q*6)*u3(i) 79 v2(i+q*0)=u2(i+q*1)*du0+u2(i+q*4)*du1+u2(i+q*7)*u3(i) 80 v2(i+q*1)=u2(i+q*2)*du0+u2(i+q*5)*du1+u2(i+q*8)*u3(i) 81 enddo 82 83 ierr=0 84 end 85!----------------------------------------------------------------------- 86 program test 87 88 include 'ceedf.h' 89 90 integer ceed,err,i,j,k 91 integer imode 92 parameter(imode=ceed_noninterlaced) 93 integer stridesu(3),stridesqd(3) 94 integer erestrictx,erestrictu,erestrictui 95 integer erestrictqi,erestrictlini 96 integer bx,bu 97 integer qf_setup_mass,qf_setup_diff,qf_apply,qf_apply_lin 98 integer op_setup_mass,op_setup_diff,op_apply,op_apply_lin 99 integer qdata_mass,qdata_diff,x,a,u,v 100 integer nelem,p,q,d 101 integer row,col,offset 102 parameter(nelem=6) 103 parameter(p=3) 104 parameter(q=4) 105 parameter(d=2) 106 integer ndofs,nqpts,nx,ny 107 parameter(nx=3) 108 parameter(ny=2) 109 parameter(ndofs=(nx*2+1)*(ny*2+1)) 110 parameter(nqpts=nelem*q*q) 111 integer indx(nelem*p*p) 112 real*8 arrx(d*ndofs),vv(ndofs) 113 real*8 total 114 integer*8 xoffset,voffset 115 116 character arg*32 117 118 external setup_mass,setup_diff,apply,apply_lin 119 120 call getarg(1,arg) 121 122 call ceedinit(trim(arg)//char(0),ceed,err) 123 124! DoF Coordinates 125 do i=0,nx*2 126 do j=0,ny*2 127 arrx(i+j*(nx*2+1)+0*ndofs+1)=1.d0*i/(2*nx) 128 arrx(i+j*(nx*2+1)+1*ndofs+1)=1.d0*j/(2*ny) 129 enddo 130 enddo 131 call ceedvectorcreate(ceed,d*ndofs,x,err) 132 xoffset=0 133 call ceedvectorsetarray(x,ceed_mem_host,ceed_use_pointer,arrx,xoffset,err) 134 135! Qdata Vector 136 call ceedvectorcreate(ceed,nqpts,qdata_mass,err) 137 call ceedvectorcreate(ceed,nqpts*d*(d+1)/2,qdata_diff,err) 138 139! Element Setup 140 do i=0,nelem-1 141 col=mod(i,nx) 142 row=i/nx 143 offset=col*(p-1)+row*(nx*2+1)*(p-1) 144 do j=0,p-1 145 do k=0,p-1 146 indx(p*(p*i+k)+j+1)=offset+k*(nx*2+1)+j 147 enddo 148 enddo 149 enddo 150 151! Restrictions 152 call ceedelemrestrictioncreate(ceed,imode,nelem,p*p,ndofs,d,& 153 & ceed_mem_host,ceed_use_pointer,indx,erestrictx,err) 154 155 call ceedelemrestrictioncreate(ceed,imode,nelem,p*p,ndofs,1,& 156 & ceed_mem_host,ceed_use_pointer,indx,erestrictu,err) 157 stridesu=[1,q*q,q*q] 158 call ceedelemrestrictioncreatestrided(ceed,nelem,q*q,nqpts,& 159 & 1,stridesu,erestrictui,err) 160 161 stridesqd=[1,q*q,q*q*d*(d+1)/2] 162 call ceedelemrestrictioncreatestrided(ceed,nelem,q*q,nqpts,& 163 & d*(d+1)/2,stridesqd,erestrictqi,err) 164 165! Bases 166 call ceedbasiscreatetensorh1lagrange(ceed,d,d,p,q,ceed_gauss,bx,err) 167 call ceedbasiscreatetensorh1lagrange(ceed,d,1,p,q,ceed_gauss,bu,err) 168 169! QFunction - setup mass 170 call ceedqfunctioncreateinterior(ceed,1,setup_mass,& 171 &SOURCE_DIR& 172 &//'t532-operator.h:setup_mass'//char(0),qf_setup_mass,err) 173 call ceedqfunctionaddinput(qf_setup_mass,'dx',d*d,ceed_eval_grad,err) 174 call ceedqfunctionaddinput(qf_setup_mass,'_weight',1,ceed_eval_weight,err) 175 call ceedqfunctionaddoutput(qf_setup_mass,'qdata',1,ceed_eval_none,err) 176 177! Operator - setup mass 178 call ceedoperatorcreate(ceed,qf_setup_mass,ceed_qfunction_none,& 179 & ceed_qfunction_none,op_setup_mass,err) 180 call ceedoperatorsetfield(op_setup_mass,'dx',erestrictx,& 181 & bx,ceed_vector_active,err) 182 call ceedoperatorsetfield(op_setup_mass,'_weight',& 183 & ceed_elemrestriction_none,bx,ceed_vector_none,err) 184 call ceedoperatorsetfield(op_setup_mass,'qdata',erestrictui,& 185 & ceed_basis_collocated,ceed_vector_active,err) 186 187! QFunction - setup diff 188 call ceedqfunctioncreateinterior(ceed,1,setup_diff,& 189 &SOURCE_DIR& 190 &//'t532-operator.h:setup_diff'//char(0),qf_setup_diff,err) 191 call ceedqfunctionaddinput(qf_setup_diff,'dx',d*d,ceed_eval_grad,err) 192 call ceedqfunctionaddinput(qf_setup_diff,'_weight',1,ceed_eval_weight,err) 193 call ceedqfunctionaddoutput(qf_setup_diff,'qdata',& 194 & d*(d+1)/2,ceed_eval_none,err) 195 196! Operator - setup diff 197 call ceedoperatorcreate(ceed,qf_setup_diff,ceed_qfunction_none,& 198 & ceed_qfunction_none,op_setup_diff,err) 199 call ceedoperatorsetfield(op_setup_diff,'dx',erestrictx,& 200 & bx,ceed_vector_active,err) 201 call ceedoperatorsetfield(op_setup_diff,'_weight',& 202 & ceed_elemrestriction_none,bx,ceed_vector_none,err) 203 call ceedoperatorsetfield(op_setup_diff,'qdata',erestrictqi,& 204 & ceed_basis_collocated,ceed_vector_active,err) 205 206! Apply Setup Operators 207 call ceedoperatorapply(op_setup_mass,x,qdata_mass,& 208 & ceed_request_immediate,err) 209 call ceedoperatorapply(op_setup_diff,x,qdata_diff,& 210 & ceed_request_immediate,err) 211 212! QFunction - apply 213 call ceedqfunctioncreateinterior(ceed,1,apply,& 214 &SOURCE_DIR& 215 &//'t532-operator.h:apply'//char(0),qf_apply,err) 216 call ceedqfunctionaddinput(qf_apply,'du',d,ceed_eval_grad,err) 217 call ceedqfunctionaddinput(qf_apply,'qdata_mass',1,ceed_eval_none,err) 218 call ceedqfunctionaddinput(qf_apply,'qdata_diff',& 219 & d*(d+1)/2,ceed_eval_none,err) 220 call ceedqfunctionaddinput(qf_apply,'u',1,ceed_eval_interp,err) 221 call ceedqfunctionaddoutput(qf_apply,'v',1,ceed_eval_interp,err) 222 call ceedqfunctionaddoutput(qf_apply,'dv',d,ceed_eval_grad,err) 223 224! Operator - apply 225 call ceedoperatorcreate(ceed,qf_apply,ceed_qfunction_none,& 226 & ceed_qfunction_none,op_apply,err) 227 call ceedoperatorsetfield(op_apply,'du',erestrictu,& 228 & bu,ceed_vector_active,err) 229 call ceedoperatorsetfield(op_apply,'qdata_mass',erestrictui,& 230 & ceed_basis_collocated,qdata_mass,err) 231 call ceedoperatorsetfield(op_apply,'qdata_diff',erestrictqi,& 232 & ceed_basis_collocated,qdata_diff,err) 233 call ceedoperatorsetfield(op_apply,'u',erestrictu,& 234 & bu,ceed_vector_active,err) 235 call ceedoperatorsetfield(op_apply,'v',erestrictu,& 236 & bu,ceed_vector_active,err) 237 call ceedoperatorsetfield(op_apply,'dv',erestrictu,& 238 & bu,ceed_vector_active,err) 239 240! Apply Original Operator 241 call ceedvectorcreate(ceed,ndofs,u,err) 242 call ceedvectorsetvalue(u,1.d0,err) 243 call ceedvectorcreate(ceed,ndofs,v,err) 244 call ceedvectorsetvalue(v,0.d0,err) 245 call ceedoperatorapply(op_apply,u,v,ceed_request_immediate,err) 246 247! Check Output 248 call ceedvectorgetarrayread(v,ceed_mem_host,vv,voffset,err) 249 total=0. 250 do i=1,ndofs 251 total=total+vv(voffset+i) 252 enddo 253 if (abs(total-1.)>1.0d-10) then 254! LCOV_EXCL_START 255 write(*,*) 'Error: True operator computed area = ',total,' != 1.0' 256! LCOV_EXCL_STOP 257 endif 258 call ceedvectorrestorearrayread(v,vv,voffset,err) 259 260! Assemble QFunction 261 call ceedoperatorassemblelinearqfunction(op_apply,a,erestrictlini,& 262 & ceed_request_immediate,err) 263 264! QFunction - apply linearized 265 call ceedqfunctioncreateinterior(ceed,1,apply_lin,& 266 &SOURCE_DIR& 267 &//'t532-operator.h:apply_lin'//char(0),qf_apply_lin,err) 268 call ceedqfunctionaddinput(qf_apply_lin,'du',d,ceed_eval_grad,err) 269 call ceedqfunctionaddinput(qf_apply_lin,'qdata',(d+1)*(d+1),& 270 & ceed_eval_none,err) 271 call ceedqfunctionaddinput(qf_apply_lin,'u',1,ceed_eval_interp,err) 272 call ceedqfunctionaddoutput(qf_apply_lin,'v',1,ceed_eval_interp,err) 273 call ceedqfunctionaddoutput(qf_apply_lin,'dv',d,ceed_eval_grad,err) 274 275! Operator - apply linearized 276 call ceedoperatorcreate(ceed,qf_apply_lin,ceed_qfunction_none,& 277 & ceed_qfunction_none,op_apply_lin,err) 278 call ceedoperatorsetfield(op_apply_lin,'du',erestrictu,& 279 & bu,ceed_vector_active,err) 280 call ceedoperatorsetfield(op_apply_lin,'qdata',erestrictlini,& 281 & ceed_basis_collocated,a,err) 282 call ceedoperatorsetfield(op_apply_lin,'u',erestrictu,& 283 & bu,ceed_vector_active,err) 284 call ceedoperatorsetfield(op_apply_lin,'v',erestrictu,& 285 & bu,ceed_vector_active,err) 286 call ceedoperatorsetfield(op_apply_lin,'dv',erestrictu,& 287 & bu,ceed_vector_active,err) 288 289! Apply Linearized QFunction Operator 290 call ceedvectorsetvalue(v,0.d0,err) 291 call ceedoperatorapply(op_apply_lin,u,v,ceed_request_immediate,err) 292 293! Check Output 294 call ceedvectorgetarrayread(v,ceed_mem_host,vv,voffset,err) 295 total=0. 296 do i=1,ndofs 297 total=total+vv(voffset+i) 298 enddo 299 if (abs(total-1.)>1.0d-10) then 300! LCOV_EXCL_START 301 write(*,*) 'Error: Assembled operator computed area = ',total,' != 1.0' 302! LCOV_EXCL_STOP 303 endif 304 call ceedvectorrestorearrayread(v,vv,voffset,err) 305 306! Cleanup 307 call ceedqfunctiondestroy(qf_setup_mass,err) 308 call ceedqfunctiondestroy(qf_setup_diff,err) 309 call ceedqfunctiondestroy(qf_apply,err) 310 call ceedqfunctiondestroy(qf_apply_lin,err) 311 call ceedoperatordestroy(op_setup_mass,err) 312 call ceedoperatordestroy(op_setup_diff,err) 313 call ceedoperatordestroy(op_apply,err) 314 call ceedoperatordestroy(op_apply_lin,err) 315 call ceedelemrestrictiondestroy(erestrictu,err) 316 call ceedelemrestrictiondestroy(erestrictx,err) 317 call ceedelemrestrictiondestroy(erestrictui,err) 318 call ceedelemrestrictiondestroy(erestrictqi,err) 319 call ceedelemrestrictiondestroy(erestrictlini,err) 320 call ceedbasisdestroy(bu,err) 321 call ceedbasisdestroy(bx,err) 322 call ceedvectordestroy(x,err) 323 call ceedvectordestroy(a,err) 324 call ceedvectordestroy(u,err) 325 call ceedvectordestroy(v,err) 326 call ceedvectordestroy(qdata_mass,err) 327 call ceedvectordestroy(qdata_diff,err) 328 call ceeddestroy(ceed,err) 329 end 330!----------------------------------------------------------------------- 331