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