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 erestrictx,erestrictu,erestrictxi,erestrictui 94 integer erestrictqi,erestrictlini 95 integer bx,bu 96 integer qf_setup_mass,qf_setup_diff,qf_apply,qf_apply_lin 97 integer op_setup_mass,op_setup_diff,op_apply,op_apply_lin 98 integer qdata_mass,qdata_diff,x,a,u,v 99 integer nelem,p,q,d 100 integer row,col,offset 101 parameter(nelem=6) 102 parameter(p=3) 103 parameter(q=4) 104 parameter(d=2) 105 integer ndofs,nqpts,nx,ny 106 parameter(nx=3) 107 parameter(ny=2) 108 parameter(ndofs=(nx*2+1)*(ny*2+1)) 109 parameter(nqpts=nelem*q*q) 110 integer indx(nelem*p*p) 111 real*8 arrx(d*ndofs),vv(ndofs) 112 real*8 total 113 integer*8 xoffset,voffset 114 115 character arg*32 116 117 external setup_mass,setup_diff,apply,apply_lin 118 119 call getarg(1,arg) 120 121 call ceedinit(trim(arg)//char(0),ceed,err) 122 123! DoF Coordinates 124 do i=0,nx*2 125 do j=0,ny*2 126 arrx(i+j*(nx*2+1)+0*ndofs+1)=1.d0*i/(2*nx) 127 arrx(i+j*(nx*2+1)+1*ndofs+1)=1.d0*j/(2*ny) 128 enddo 129 enddo 130 call ceedvectorcreate(ceed,d*ndofs,x,err) 131 xoffset=0 132 call ceedvectorsetarray(x,ceed_mem_host,ceed_use_pointer,arrx,xoffset,err) 133 134! Qdata Vector 135 call ceedvectorcreate(ceed,nqpts,qdata_mass,err) 136 call ceedvectorcreate(ceed,nqpts*d*(d+1)/2,qdata_diff,err) 137 138! Element Setup 139 do i=0,nelem-1 140 col=mod(i,nx) 141 row=i/nx 142 offset=col*(p-1)+row*(nx*2+1)*(p-1) 143 do j=0,p-1 144 do k=0,p-1 145 indx(p*(p*i+k)+j+1)=offset+k*(nx*2+1)+j 146 enddo 147 enddo 148 enddo 149 150! Restrictions 151 call ceedelemrestrictioncreate(ceed,imode,nelem,p*p,ndofs,d,& 152 & ceed_mem_host,ceed_use_pointer,indx,erestrictx,err) 153 call ceedelemrestrictioncreateidentity(ceed,imode,nelem,p*p,& 154 & nelem*p*p,d,erestrictxi,err) 155 156 call ceedelemrestrictioncreate(ceed,imode,nelem,p*p,ndofs,1,& 157 & ceed_mem_host,ceed_use_pointer,indx,erestrictu,err) 158 call ceedelemrestrictioncreateidentity(ceed,imode,nelem,q*q,nqpts,& 159 & 1,erestrictui,err) 160 161 call ceedelemrestrictioncreateidentity(ceed,imode,nelem,q*q,nqpts,& 162 & d*(d+1)/2,erestrictqi,err) 163 164! Bases 165 call ceedbasiscreatetensorh1lagrange(ceed,d,d,p,q,ceed_gauss,bx,err) 166 call ceedbasiscreatetensorh1lagrange(ceed,d,1,p,q,ceed_gauss,bu,err) 167 168! QFunction - setup mass 169 call ceedqfunctioncreateinterior(ceed,1,setup_mass,& 170 &SOURCE_DIR& 171 &//'t532-operator.h:setup_mass'//char(0),qf_setup_mass,err) 172 call ceedqfunctionaddinput(qf_setup_mass,'dx',d*d,ceed_eval_grad,err) 173 call ceedqfunctionaddinput(qf_setup_mass,'_weight',1,ceed_eval_weight,err) 174 call ceedqfunctionaddoutput(qf_setup_mass,'qdata',1,ceed_eval_none,err) 175 176! Operator - setup mass 177 call ceedoperatorcreate(ceed,qf_setup_mass,ceed_qfunction_none,& 178 & ceed_qfunction_none,op_setup_mass,err) 179 call ceedoperatorsetfield(op_setup_mass,'dx',erestrictx,& 180 & bx,ceed_vector_active,err) 181 call ceedoperatorsetfield(op_setup_mass,'_weight',erestrictxi,& 182 & bx,ceed_vector_none,err) 183 call ceedoperatorsetfield(op_setup_mass,'qdata',erestrictui,& 184 & ceed_basis_collocated,ceed_vector_active,err) 185 186! QFunction - setup diff 187 call ceedqfunctioncreateinterior(ceed,1,setup_diff,& 188 &SOURCE_DIR& 189 &//'t532-operator.h:setup_diff'//char(0),qf_setup_diff,err) 190 call ceedqfunctionaddinput(qf_setup_diff,'dx',d*d,ceed_eval_grad,err) 191 call ceedqfunctionaddinput(qf_setup_diff,'_weight',1,ceed_eval_weight,err) 192 call ceedqfunctionaddoutput(qf_setup_diff,'qdata',& 193 & d*(d+1)/2,ceed_eval_none,err) 194 195! Operator - setup diff 196 call ceedoperatorcreate(ceed,qf_setup_diff,ceed_qfunction_none,& 197 & ceed_qfunction_none,op_setup_diff,err) 198 call ceedoperatorsetfield(op_setup_diff,'dx',erestrictx,& 199 & bx,ceed_vector_active,err) 200 call ceedoperatorsetfield(op_setup_diff,'_weight',erestrictxi,& 201 & bx,ceed_vector_none,err) 202 call ceedoperatorsetfield(op_setup_diff,'qdata',erestrictqi,& 203 & ceed_basis_collocated,ceed_vector_active,err) 204 205! Apply Setup Operators 206 call ceedoperatorapply(op_setup_mass,x,qdata_mass,& 207 & ceed_request_immediate,err) 208 call ceedoperatorapply(op_setup_diff,x,qdata_diff,& 209 & ceed_request_immediate,err) 210 211! QFunction - apply 212 call ceedqfunctioncreateinterior(ceed,1,apply,& 213 &SOURCE_DIR& 214 &//'t532-operator.h:apply'//char(0),qf_apply,err) 215 call ceedqfunctionaddinput(qf_apply,'du',d,ceed_eval_grad,err) 216 call ceedqfunctionaddinput(qf_apply,'qdata_mass',1,ceed_eval_none,err) 217 call ceedqfunctionaddinput(qf_apply,'qdata_diff',& 218 & d*(d+1)/2,ceed_eval_none,err) 219 call ceedqfunctionaddinput(qf_apply,'u',1,ceed_eval_interp,err) 220 call ceedqfunctionaddoutput(qf_apply,'v',1,ceed_eval_interp,err) 221 call ceedqfunctionaddoutput(qf_apply,'dv',d,ceed_eval_grad,err) 222 223! Operator - apply 224 call ceedoperatorcreate(ceed,qf_apply,ceed_qfunction_none,& 225 & ceed_qfunction_none,op_apply,err) 226 call ceedoperatorsetfield(op_apply,'du',erestrictu,& 227 & bu,ceed_vector_active,err) 228 call ceedoperatorsetfield(op_apply,'qdata_mass',erestrictui,& 229 & ceed_basis_collocated,qdata_mass,err) 230 call ceedoperatorsetfield(op_apply,'qdata_diff',erestrictqi,& 231 & ceed_basis_collocated,qdata_diff,err) 232 call ceedoperatorsetfield(op_apply,'u',erestrictu,& 233 & bu,ceed_vector_active,err) 234 call ceedoperatorsetfield(op_apply,'v',erestrictu,& 235 & bu,ceed_vector_active,err) 236 call ceedoperatorsetfield(op_apply,'dv',erestrictu,& 237 & bu,ceed_vector_active,err) 238 239! Apply Original Operator 240 call ceedvectorcreate(ceed,ndofs,u,err) 241 call ceedvectorsetvalue(u,1.d0,err) 242 call ceedvectorcreate(ceed,ndofs,v,err) 243 call ceedvectorsetvalue(v,0.d0,err) 244 call ceedoperatorapply(op_apply,u,v,ceed_request_immediate,err) 245 246! Check Output 247 call ceedvectorgetarrayread(v,ceed_mem_host,vv,voffset,err) 248 total=0. 249 do i=1,ndofs 250 total=total+vv(voffset+i) 251 enddo 252 if (abs(total-1.)>1.0d-10) then 253! LCOV_EXCL_START 254 write(*,*) 'Error: True operator computed area = ',total,' != 1.0' 255! LCOV_EXCL_STOP 256 endif 257 call ceedvectorrestorearrayread(v,vv,voffset,err) 258 259! Assemble QFunction 260 call ceedoperatorassemblelinearqfunction(op_apply,a,erestrictlini,& 261 & ceed_request_immediate,err) 262 263! QFunction - apply linearized 264 call ceedqfunctioncreateinterior(ceed,1,apply_lin,& 265 &SOURCE_DIR& 266 &//'t532-operator.h:apply_lin'//char(0),qf_apply_lin,err) 267 call ceedqfunctionaddinput(qf_apply_lin,'du',d,ceed_eval_grad,err) 268 call ceedqfunctionaddinput(qf_apply_lin,'qdata',(d+1)*(d+1),& 269 & ceed_eval_none,err) 270 call ceedqfunctionaddinput(qf_apply_lin,'u',1,ceed_eval_interp,err) 271 call ceedqfunctionaddoutput(qf_apply_lin,'v',1,ceed_eval_interp,err) 272 call ceedqfunctionaddoutput(qf_apply_lin,'dv',d,ceed_eval_grad,err) 273 274! Operator - apply linearized 275 call ceedoperatorcreate(ceed,qf_apply_lin,ceed_qfunction_none,& 276 & ceed_qfunction_none,op_apply_lin,err) 277 call ceedoperatorsetfield(op_apply_lin,'du',erestrictu,& 278 & bu,ceed_vector_active,err) 279 call ceedoperatorsetfield(op_apply_lin,'qdata',erestrictlini,& 280 & ceed_basis_collocated,a,err) 281 call ceedoperatorsetfield(op_apply_lin,'u',erestrictu,& 282 & bu,ceed_vector_active,err) 283 call ceedoperatorsetfield(op_apply_lin,'v',erestrictu,& 284 & bu,ceed_vector_active,err) 285 call ceedoperatorsetfield(op_apply_lin,'dv',erestrictu,& 286 & bu,ceed_vector_active,err) 287 288! Apply Linearized QFunction Operator 289 call ceedvectorsetvalue(v,0.d0,err) 290 call ceedoperatorapply(op_apply_lin,u,v,ceed_request_immediate,err) 291 292! Check Output 293 call ceedvectorgetarrayread(v,ceed_mem_host,vv,voffset,err) 294 total=0. 295 do i=1,ndofs 296 total=total+vv(voffset+i) 297 enddo 298 if (abs(total-1.)>1.0d-10) then 299! LCOV_EXCL_START 300 write(*,*) 'Error: Assembled operator computed area = ',total,' != 1.0' 301! LCOV_EXCL_STOP 302 endif 303 call ceedvectorrestorearrayread(v,vv,voffset,err) 304 305! Cleanup 306 call ceedqfunctiondestroy(qf_setup_mass,err) 307 call ceedqfunctiondestroy(qf_setup_diff,err) 308 call ceedqfunctiondestroy(qf_apply,err) 309 call ceedqfunctiondestroy(qf_apply_lin,err) 310 call ceedoperatordestroy(op_setup_mass,err) 311 call ceedoperatordestroy(op_setup_diff,err) 312 call ceedoperatordestroy(op_apply,err) 313 call ceedoperatordestroy(op_apply_lin,err) 314 call ceedelemrestrictiondestroy(erestrictu,err) 315 call ceedelemrestrictiondestroy(erestrictx,err) 316 call ceedelemrestrictiondestroy(erestrictxi,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