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 stridesx(3),stridesu(3),stridesqd(3) 94 integer erestrictx,erestrictu,erestrictxi,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 stridesx=[1,p*p,p*p*d] 155 call ceedelemrestrictioncreatestrided(ceed,nelem,p*p,& 156 & nelem*p*p,d,stridesx,erestrictxi,err) 157 158 call ceedelemrestrictioncreate(ceed,imode,nelem,p*p,ndofs,1,& 159 & ceed_mem_host,ceed_use_pointer,indx,erestrictu,err) 160 stridesu=[1,q*q,q*q] 161 call ceedelemrestrictioncreatestrided(ceed,nelem,q*q,nqpts,& 162 & 1,stridesu,erestrictui,err) 163 164 stridesqd=[1,q*q,q*q*d*(d+1)/2] 165 call ceedelemrestrictioncreatestrided(ceed,nelem,q*q,nqpts,& 166 & d*(d+1)/2,stridesqd,erestrictqi,err) 167 168! Bases 169 call ceedbasiscreatetensorh1lagrange(ceed,d,d,p,q,ceed_gauss,bx,err) 170 call ceedbasiscreatetensorh1lagrange(ceed,d,1,p,q,ceed_gauss,bu,err) 171 172! QFunction - setup mass 173 call ceedqfunctioncreateinterior(ceed,1,setup_mass,& 174 &SOURCE_DIR& 175 &//'t532-operator.h:setup_mass'//char(0),qf_setup_mass,err) 176 call ceedqfunctionaddinput(qf_setup_mass,'dx',d*d,ceed_eval_grad,err) 177 call ceedqfunctionaddinput(qf_setup_mass,'_weight',1,ceed_eval_weight,err) 178 call ceedqfunctionaddoutput(qf_setup_mass,'qdata',1,ceed_eval_none,err) 179 180! Operator - setup mass 181 call ceedoperatorcreate(ceed,qf_setup_mass,ceed_qfunction_none,& 182 & ceed_qfunction_none,op_setup_mass,err) 183 call ceedoperatorsetfield(op_setup_mass,'dx',erestrictx,& 184 & bx,ceed_vector_active,err) 185 call ceedoperatorsetfield(op_setup_mass,'_weight',erestrictxi,& 186 & bx,ceed_vector_none,err) 187 call ceedoperatorsetfield(op_setup_mass,'qdata',erestrictui,& 188 & ceed_basis_collocated,ceed_vector_active,err) 189 190! QFunction - setup diff 191 call ceedqfunctioncreateinterior(ceed,1,setup_diff,& 192 &SOURCE_DIR& 193 &//'t532-operator.h:setup_diff'//char(0),qf_setup_diff,err) 194 call ceedqfunctionaddinput(qf_setup_diff,'dx',d*d,ceed_eval_grad,err) 195 call ceedqfunctionaddinput(qf_setup_diff,'_weight',1,ceed_eval_weight,err) 196 call ceedqfunctionaddoutput(qf_setup_diff,'qdata',& 197 & d*(d+1)/2,ceed_eval_none,err) 198 199! Operator - setup diff 200 call ceedoperatorcreate(ceed,qf_setup_diff,ceed_qfunction_none,& 201 & ceed_qfunction_none,op_setup_diff,err) 202 call ceedoperatorsetfield(op_setup_diff,'dx',erestrictx,& 203 & bx,ceed_vector_active,err) 204 call ceedoperatorsetfield(op_setup_diff,'_weight',erestrictxi,& 205 & bx,ceed_vector_none,err) 206 call ceedoperatorsetfield(op_setup_diff,'qdata',erestrictqi,& 207 & ceed_basis_collocated,ceed_vector_active,err) 208 209! Apply Setup Operators 210 call ceedoperatorapply(op_setup_mass,x,qdata_mass,& 211 & ceed_request_immediate,err) 212 call ceedoperatorapply(op_setup_diff,x,qdata_diff,& 213 & ceed_request_immediate,err) 214 215! QFunction - apply 216 call ceedqfunctioncreateinterior(ceed,1,apply,& 217 &SOURCE_DIR& 218 &//'t532-operator.h:apply'//char(0),qf_apply,err) 219 call ceedqfunctionaddinput(qf_apply,'du',d,ceed_eval_grad,err) 220 call ceedqfunctionaddinput(qf_apply,'qdata_mass',1,ceed_eval_none,err) 221 call ceedqfunctionaddinput(qf_apply,'qdata_diff',& 222 & d*(d+1)/2,ceed_eval_none,err) 223 call ceedqfunctionaddinput(qf_apply,'u',1,ceed_eval_interp,err) 224 call ceedqfunctionaddoutput(qf_apply,'v',1,ceed_eval_interp,err) 225 call ceedqfunctionaddoutput(qf_apply,'dv',d,ceed_eval_grad,err) 226 227! Operator - apply 228 call ceedoperatorcreate(ceed,qf_apply,ceed_qfunction_none,& 229 & ceed_qfunction_none,op_apply,err) 230 call ceedoperatorsetfield(op_apply,'du',erestrictu,& 231 & bu,ceed_vector_active,err) 232 call ceedoperatorsetfield(op_apply,'qdata_mass',erestrictui,& 233 & ceed_basis_collocated,qdata_mass,err) 234 call ceedoperatorsetfield(op_apply,'qdata_diff',erestrictqi,& 235 & ceed_basis_collocated,qdata_diff,err) 236 call ceedoperatorsetfield(op_apply,'u',erestrictu,& 237 & bu,ceed_vector_active,err) 238 call ceedoperatorsetfield(op_apply,'v',erestrictu,& 239 & bu,ceed_vector_active,err) 240 call ceedoperatorsetfield(op_apply,'dv',erestrictu,& 241 & bu,ceed_vector_active,err) 242 243! Apply Original Operator 244 call ceedvectorcreate(ceed,ndofs,u,err) 245 call ceedvectorsetvalue(u,1.d0,err) 246 call ceedvectorcreate(ceed,ndofs,v,err) 247 call ceedvectorsetvalue(v,0.d0,err) 248 call ceedoperatorapply(op_apply,u,v,ceed_request_immediate,err) 249 250! Check Output 251 call ceedvectorgetarrayread(v,ceed_mem_host,vv,voffset,err) 252 total=0. 253 do i=1,ndofs 254 total=total+vv(voffset+i) 255 enddo 256 if (abs(total-1.)>1.0d-10) then 257! LCOV_EXCL_START 258 write(*,*) 'Error: True operator computed area = ',total,' != 1.0' 259! LCOV_EXCL_STOP 260 endif 261 call ceedvectorrestorearrayread(v,vv,voffset,err) 262 263! Assemble QFunction 264 call ceedoperatorassemblelinearqfunction(op_apply,a,erestrictlini,& 265 & ceed_request_immediate,err) 266 267! QFunction - apply linearized 268 call ceedqfunctioncreateinterior(ceed,1,apply_lin,& 269 &SOURCE_DIR& 270 &//'t532-operator.h:apply_lin'//char(0),qf_apply_lin,err) 271 call ceedqfunctionaddinput(qf_apply_lin,'du',d,ceed_eval_grad,err) 272 call ceedqfunctionaddinput(qf_apply_lin,'qdata',(d+1)*(d+1),& 273 & ceed_eval_none,err) 274 call ceedqfunctionaddinput(qf_apply_lin,'u',1,ceed_eval_interp,err) 275 call ceedqfunctionaddoutput(qf_apply_lin,'v',1,ceed_eval_interp,err) 276 call ceedqfunctionaddoutput(qf_apply_lin,'dv',d,ceed_eval_grad,err) 277 278! Operator - apply linearized 279 call ceedoperatorcreate(ceed,qf_apply_lin,ceed_qfunction_none,& 280 & ceed_qfunction_none,op_apply_lin,err) 281 call ceedoperatorsetfield(op_apply_lin,'du',erestrictu,& 282 & bu,ceed_vector_active,err) 283 call ceedoperatorsetfield(op_apply_lin,'qdata',erestrictlini,& 284 & ceed_basis_collocated,a,err) 285 call ceedoperatorsetfield(op_apply_lin,'u',erestrictu,& 286 & bu,ceed_vector_active,err) 287 call ceedoperatorsetfield(op_apply_lin,'v',erestrictu,& 288 & bu,ceed_vector_active,err) 289 call ceedoperatorsetfield(op_apply_lin,'dv',erestrictu,& 290 & bu,ceed_vector_active,err) 291 292! Apply Linearized QFunction Operator 293 call ceedvectorsetvalue(v,0.d0,err) 294 call ceedoperatorapply(op_apply_lin,u,v,ceed_request_immediate,err) 295 296! Check Output 297 call ceedvectorgetarrayread(v,ceed_mem_host,vv,voffset,err) 298 total=0. 299 do i=1,ndofs 300 total=total+vv(voffset+i) 301 enddo 302 if (abs(total-1.)>1.0d-10) then 303! LCOV_EXCL_START 304 write(*,*) 'Error: Assembled operator computed area = ',total,' != 1.0' 305! LCOV_EXCL_STOP 306 endif 307 call ceedvectorrestorearrayread(v,vv,voffset,err) 308 309! Cleanup 310 call ceedqfunctiondestroy(qf_setup_mass,err) 311 call ceedqfunctiondestroy(qf_setup_diff,err) 312 call ceedqfunctiondestroy(qf_apply,err) 313 call ceedqfunctiondestroy(qf_apply_lin,err) 314 call ceedoperatordestroy(op_setup_mass,err) 315 call ceedoperatordestroy(op_setup_diff,err) 316 call ceedoperatordestroy(op_apply,err) 317 call ceedoperatordestroy(op_apply_lin,err) 318 call ceedelemrestrictiondestroy(erestrictu,err) 319 call ceedelemrestrictiondestroy(erestrictx,err) 320 call ceedelemrestrictiondestroy(erestrictxi,err) 321 call ceedelemrestrictiondestroy(erestrictui,err) 322 call ceedelemrestrictiondestroy(erestrictqi,err) 323 call ceedelemrestrictiondestroy(erestrictlini,err) 324 call ceedbasisdestroy(bu,err) 325 call ceedbasisdestroy(bx,err) 326 call ceedvectordestroy(x,err) 327 call ceedvectordestroy(a,err) 328 call ceedvectordestroy(u,err) 329 call ceedvectordestroy(v,err) 330 call ceedvectordestroy(qdata_mass,err) 331 call ceedvectordestroy(qdata_diff,err) 332 call ceeddestroy(ceed,err) 333 end 334!----------------------------------------------------------------------- 335