1!----------------------------------------------------------------------- 2 subroutine setup(ctx,q,u1,u2,u3,u4,u5,u6,u7,u8,u9,u10,u11,u12,u13,u14,& 3& u15,u16,v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,v13,v14,v15,v16,ierr) 4 real*8 ctx 5 real*8 u1(1) 6 real*8 u2(1) 7 real*8 v1(1) 8 real*8 w 9 integer q,ierr 10 11 do i=1,q 12 w=u2(i)/(u1(i+q*0)*u1(i+q*3)-u1(i+q*1)*u1(i+q*2)) 13 v1(i+q*0)=w*(u1(i+q*2)*u1(i+q*2)+u1(i+q*3)*u1(i+q*3)) 14 v1(i+q*1)=-w*(u1(i+q*0)*u1(i+q*2)+u1(i+q*2)*u1(i+q*3)) 15 v1(i+q*2)=w*(u1(i+q*0)*u1(i+q*0)+u1(i+q*1)*u1(i+q*1)) 16 enddo 17 18 ierr=0 19 end 20!----------------------------------------------------------------------- 21 subroutine diff(ctx,q,u1,u2,u3,u4,u5,u6,u7,u8,u9,u10,u11,u12,u13,u14,& 22& u15,u16,v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,v13,v14,v15,v16,ierr) 23 real*8 ctx 24 real*8 u1(1) 25 real*8 u2(1) 26 real*8 v1(1) 27 real*8 du0,du1 28 integer q,ierr 29 30 do i=1,q 31 du0=u1(i+q*0) 32 du1=u1(i+q*1) 33 v1(i+q*0)=u2(i+q*0)*du0+u2(i+q*1)*du1 34 v1(i+q*1)=u2(i+q*1)*du0+u2(i+q*2)*du1 35 enddo 36 37 ierr=0 38 end 39!----------------------------------------------------------------------- 40 subroutine diff_lin(ctx,q,u1,u2,u3,u4,u5,u6,u7,u8,u9,u10,u11,u12,u13,u14,& 41& u15,u16,v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,v13,v14,v15,v16,ierr) 42 real*8 ctx 43 real*8 u1(1) 44 real*8 u2(1) 45 real*8 v1(1) 46 real*8 du0,du1 47 integer q,ierr 48 49 do i=1,q 50 du0=u1(i+q*0) 51 du1=u1(i+q*1) 52 v1(i+q*0)=u2(i+q*0)*du0+u2(i+q*1)*du1 53 v1(i+q*1)=u2(i+q*2)*du0+u2(i+q*3)*du1 54 enddo 55 56 ierr=0 57 end 58!----------------------------------------------------------------------- 59 program test 60 61 include 'ceedf.h' 62 63 integer ceed,err,i,j,k 64 integer erestrictx,erestrictu,erestrictxi,erestrictui 65 integer erestrictqi,erestrictlini 66 integer bx,bu 67 integer qf_setup,qf_diff,qf_diff_lin 68 integer op_setup,op_diff,op_diff_lin 69 integer qdata,x,a,u,v 70 integer nelem,p,q,d 71 integer row,col,offset 72 parameter(nelem=6) 73 parameter(p=3) 74 parameter(q=4) 75 parameter(d=2) 76 integer ndofs,nqpts,nx,ny 77 parameter(nx=3) 78 parameter(ny=2) 79 parameter(ndofs=(nx*2+1)*(ny*2+1)) 80 parameter(nqpts=nelem*q*q) 81 integer indx(nelem*p*p) 82 real*8 arrx(d*ndofs),vv(ndofs) 83 integer*8 xoffset,voffset 84 85 character arg*32 86 87 external setup,diff,diff_lin 88 89 call getarg(1,arg) 90 91 call ceedinit(trim(arg)//char(0),ceed,err) 92 93! DoF Coordinates 94 do i=0,nx*2 95 do j=0,ny*2 96 arrx(i+j*(nx*2+1)+0*ndofs+1)=1.d0*i/(2*nx) 97 arrx(i+j*(nx*2+1)+1*ndofs+1)=1.d0*j/(2*ny) 98 enddo 99 enddo 100 call ceedvectorcreate(ceed,d*ndofs,x,err) 101 xoffset=0 102 call ceedvectorsetarray(x,ceed_mem_host,ceed_use_pointer,arrx,xoffset,err) 103 104! Qdata Vector 105 call ceedvectorcreate(ceed,nqpts*d*(d+1)/2,qdata,err) 106 107! Element Setup 108 do i=0,nelem-1 109 col=mod(i,nx) 110 row=i/nx 111 offset=col*(p-1)+row*(nx*2+1)*(p-1) 112 do j=0,p-1 113 do k=0,p-1 114 indx(p*(p*i+k)+j+1)=offset+k*(nx*2+1)+j 115 enddo 116 enddo 117 enddo 118 119! Restrictions 120 call ceedelemrestrictioncreate(ceed,nelem,p*p,ndofs,d,& 121 & ceed_mem_host,ceed_use_pointer,indx,erestrictx,err) 122 call ceedelemrestrictioncreateidentity(ceed,nelem,p*p,& 123 & nelem*p*p,d,erestrictxi,err) 124 125 call ceedelemrestrictioncreate(ceed,nelem,p*p,ndofs,1,& 126 & ceed_mem_host,ceed_use_pointer,indx,erestrictu,err) 127 call ceedelemrestrictioncreateidentity(ceed,nelem,q*q,nqpts,& 128 & 1,erestrictui,err) 129 130 call ceedelemrestrictioncreateidentity(ceed,nelem,q*q,nqpts,& 131 & d*(d+1)/2,erestrictqi,err) 132 133! Bases 134 call ceedbasiscreatetensorh1lagrange(ceed,d,d,p,q,ceed_gauss,& 135 & bx,err) 136 call ceedbasiscreatetensorh1lagrange(ceed,d,1,p,q,ceed_gauss,& 137 & bu,err) 138 139! QFunction - setup 140 call ceedqfunctioncreateinterior(ceed,1,setup,& 141 &SOURCE_DIR& 142 &//'t531-operator.h:setup'//char(0),qf_setup,err) 143 call ceedqfunctionaddinput(qf_setup,'dx',d*d,ceed_eval_grad,err) 144 call ceedqfunctionaddinput(qf_setup,'_weight',1,ceed_eval_weight,err) 145 call ceedqfunctionaddoutput(qf_setup,'qdata',d*(d+1)/2,ceed_eval_none,err) 146 147! Operator - setup 148 call ceedoperatorcreate(ceed,qf_setup,ceed_qfunction_none,& 149 & ceed_qfunction_none,op_setup,err) 150 call ceedoperatorsetfield(op_setup,'dx',erestrictx,& 151 & ceed_notranspose,bx,ceed_vector_active,err) 152 call ceedoperatorsetfield(op_setup,'_weight',erestrictxi,& 153 & ceed_notranspose,bx,ceed_vector_none,err) 154 call ceedoperatorsetfield(op_setup,'qdata',erestrictqi,& 155 & ceed_notranspose,ceed_basis_collocated,ceed_vector_active,err) 156 157! Apply Setup Operator 158 call ceedoperatorapply(op_setup,x,qdata,ceed_request_immediate,err) 159 160! QFunction - apply 161 call ceedqfunctioncreateinterior(ceed,1,diff,& 162 &SOURCE_DIR& 163 &//'t531-operator.h:diff'//char(0),qf_diff,err) 164 call ceedqfunctionaddinput(qf_diff,'du',d,ceed_eval_grad,err) 165 call ceedqfunctionaddinput(qf_diff,'qdata',d*(d+1)/2,ceed_eval_none,err) 166 call ceedqfunctionaddoutput(qf_diff,'dv',d,ceed_eval_grad,err) 167 168! Operator - apply 169 call ceedoperatorcreate(ceed,qf_diff,ceed_qfunction_none,& 170 & ceed_qfunction_none,op_diff,err) 171 call ceedoperatorsetfield(op_diff,'du',erestrictu,& 172 & ceed_notranspose,bu,ceed_vector_active,err) 173 call ceedoperatorsetfield(op_diff,'qdata',erestrictqi,& 174 & ceed_notranspose,ceed_basis_collocated,qdata,err) 175 call ceedoperatorsetfield(op_diff,'dv',erestrictu,& 176 & ceed_notranspose,bu,ceed_vector_active,err) 177 178! Apply original Poisson Operator 179 call ceedvectorcreate(ceed,ndofs,u,err) 180 call ceedvectorsetvalue(u,1.d0,err) 181 call ceedvectorcreate(ceed,ndofs,v,err) 182 call ceedvectorsetvalue(v,0.d0,err) 183 call ceedoperatorapply(op_diff,u,v,ceed_request_immediate,err) 184 185! Check Output 186 call ceedvectorgetarrayread(v,ceed_mem_host,vv,voffset,err) 187 do i=1,ndofs 188 if (abs(vv(voffset+i))>1.0d-14) then 189! LCOV_EXCL_START 190 write(*,*) 'Error: Operator computed v[i] = ',vv(voffset+i),' != 0.0' 191! LCOV_EXCL_STOP 192 endif 193 enddo 194 call ceedvectorrestorearrayread(v,vv,voffset,err) 195 196! Assemble QFunction 197 call ceedoperatorassemblelinearqfunction(op_diff,a,erestrictlini,& 198 & ceed_request_immediate,err) 199 200! QFunction - apply linearized 201 call ceedqfunctioncreateinterior(ceed,1,diff_lin,& 202 &SOURCE_DIR& 203 &//'t531-operator.h:diff_lin'//char(0),qf_diff_lin,err) 204 call ceedqfunctionaddinput(qf_diff_lin,'du',d,ceed_eval_grad,err) 205 call ceedqfunctionaddinput(qf_diff_lin,'qdata',d*d,ceed_eval_none,err) 206 call ceedqfunctionaddoutput(qf_diff_lin,'dv',d,ceed_eval_grad,err) 207 208! Operator - apply linearized 209 call ceedoperatorcreate(ceed,qf_diff_lin,ceed_qfunction_none,& 210 & ceed_qfunction_none,op_diff_lin,err) 211 call ceedoperatorsetfield(op_diff_lin,'du',erestrictu,& 212 & ceed_notranspose,bu,ceed_vector_active,err) 213 call ceedoperatorsetfield(op_diff_lin,'qdata',erestrictlini,& 214 & ceed_notranspose,ceed_basis_collocated,a,err) 215 call ceedoperatorsetfield(op_diff_lin,'dv',erestrictu,& 216 & ceed_notranspose,bu,ceed_vector_active,err) 217 218! Apply linearized Poisson Operator 219 call ceedvectorsetvalue(v,0.d0,err) 220 call ceedoperatorapply(op_diff_lin,u,v,ceed_request_immediate,err) 221 222! Check Output 223 call ceedvectorgetarrayread(v,ceed_mem_host,vv,voffset,err) 224 do i=1,ndofs 225 if (abs(vv(voffset+i))>1.0d-14) then 226! LCOV_EXCL_START 227 write(*,*) 'Error: Linearized operator computed v[i] = ',vv(voffset+i),& 228 & ' != 0.0' 229! LCOV_EXCL_STOP 230 endif 231 enddo 232 call ceedvectorrestorearrayread(v,vv,voffset,err) 233 234! Cleanup 235 call ceedqfunctiondestroy(qf_setup,err) 236 call ceedqfunctiondestroy(qf_diff,err) 237 call ceedqfunctiondestroy(qf_diff_lin,err) 238 call ceedoperatordestroy(op_setup,err) 239 call ceedoperatordestroy(op_diff,err) 240 call ceedoperatordestroy(op_diff_lin,err) 241 call ceedelemrestrictiondestroy(erestrictu,err) 242 call ceedelemrestrictiondestroy(erestrictx,err) 243 call ceedelemrestrictiondestroy(erestrictxi,err) 244 call ceedelemrestrictiondestroy(erestrictui,err) 245 call ceedelemrestrictiondestroy(erestrictqi,err) 246 call ceedelemrestrictiondestroy(erestrictlini,err) 247 call ceedbasisdestroy(bu,err) 248 call ceedbasisdestroy(bx,err) 249 call ceedvectordestroy(x,err) 250 call ceedvectordestroy(a,err) 251 call ceedvectordestroy(u,err) 252 call ceedvectordestroy(v,err) 253 call ceedvectordestroy(qdata,err) 254 call ceeddestroy(ceed,err) 255 end 256!----------------------------------------------------------------------- 257