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