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 imode 65 parameter(imode=ceed_noninterlaced) 66 integer erestrictx,erestrictu,erestrictxi,erestrictui 67 integer erestrictqi,erestrictlini 68 integer bx,bu 69 integer qf_setup,qf_diff,qf_diff_lin 70 integer op_setup,op_diff,op_diff_lin 71 integer qdata,x,a,u,v 72 integer nelem,p,q,d 73 integer row,col,offset 74 parameter(nelem=6) 75 parameter(p=3) 76 parameter(q=4) 77 parameter(d=2) 78 integer ndofs,nqpts,nx,ny 79 parameter(nx=3) 80 parameter(ny=2) 81 parameter(ndofs=(nx*2+1)*(ny*2+1)) 82 parameter(nqpts=nelem*q*q) 83 integer indx(nelem*p*p) 84 real*8 arrx(d*ndofs),vv(ndofs) 85 integer*8 xoffset,voffset 86 87 character arg*32 88 89 external setup,diff,diff_lin 90 91 call getarg(1,arg) 92 93 call ceedinit(trim(arg)//char(0),ceed,err) 94 95! DoF Coordinates 96 do i=0,nx*2 97 do j=0,ny*2 98 arrx(i+j*(nx*2+1)+0*ndofs+1)=1.d0*i/(2*nx) 99 arrx(i+j*(nx*2+1)+1*ndofs+1)=1.d0*j/(2*ny) 100 enddo 101 enddo 102 call ceedvectorcreate(ceed,d*ndofs,x,err) 103 xoffset=0 104 call ceedvectorsetarray(x,ceed_mem_host,ceed_use_pointer,arrx,xoffset,err) 105 106! Qdata Vector 107 call ceedvectorcreate(ceed,nqpts*d*(d+1)/2,qdata,err) 108 109! Element Setup 110 do i=0,nelem-1 111 col=mod(i,nx) 112 row=i/nx 113 offset=col*(p-1)+row*(nx*2+1)*(p-1) 114 do j=0,p-1 115 do k=0,p-1 116 indx(p*(p*i+k)+j+1)=offset+k*(nx*2+1)+j 117 enddo 118 enddo 119 enddo 120 121! Restrictions 122 call ceedelemrestrictioncreate(ceed,imode,nelem,p*p,ndofs,d,& 123 & ceed_mem_host,ceed_use_pointer,indx,erestrictx,err) 124 call ceedelemrestrictioncreateidentity(ceed,imode,nelem,p*p,& 125 & nelem*p*p,d,erestrictxi,err) 126 127 call ceedelemrestrictioncreate(ceed,imode,nelem,p*p,ndofs,1,& 128 & ceed_mem_host,ceed_use_pointer,indx,erestrictu,err) 129 call ceedelemrestrictioncreateidentity(ceed,imode,nelem,q*q,nqpts,& 130 & 1,erestrictui,err) 131 132 call ceedelemrestrictioncreateidentity(ceed,imode,nelem,q*q,nqpts,& 133 & d*(d+1)/2,erestrictqi,err) 134 135! Bases 136 call ceedbasiscreatetensorh1lagrange(ceed,d,d,p,q,ceed_gauss,& 137 & bx,err) 138 call ceedbasiscreatetensorh1lagrange(ceed,d,1,p,q,ceed_gauss,& 139 & bu,err) 140 141! QFunction - setup 142 call ceedqfunctioncreateinterior(ceed,1,setup,& 143 &SOURCE_DIR& 144 &//'t531-operator.h:setup'//char(0),qf_setup,err) 145 call ceedqfunctionaddinput(qf_setup,'dx',d*d,ceed_eval_grad,err) 146 call ceedqfunctionaddinput(qf_setup,'_weight',1,ceed_eval_weight,err) 147 call ceedqfunctionaddoutput(qf_setup,'qdata',d*(d+1)/2,ceed_eval_none,err) 148 149! Operator - setup 150 call ceedoperatorcreate(ceed,qf_setup,ceed_qfunction_none,& 151 & ceed_qfunction_none,op_setup,err) 152 call ceedoperatorsetfield(op_setup,'dx',erestrictx,& 153 & bx,ceed_vector_active,err) 154 call ceedoperatorsetfield(op_setup,'_weight',erestrictxi,& 155 & bx,ceed_vector_none,err) 156 call ceedoperatorsetfield(op_setup,'qdata',erestrictqi,& 157 & ceed_basis_collocated,ceed_vector_active,err) 158 159! Apply Setup Operator 160 call ceedoperatorapply(op_setup,x,qdata,ceed_request_immediate,err) 161 162! QFunction - apply 163 call ceedqfunctioncreateinterior(ceed,1,diff,& 164 &SOURCE_DIR& 165 &//'t531-operator.h:diff'//char(0),qf_diff,err) 166 call ceedqfunctionaddinput(qf_diff,'du',d,ceed_eval_grad,err) 167 call ceedqfunctionaddinput(qf_diff,'qdata',d*(d+1)/2,ceed_eval_none,err) 168 call ceedqfunctionaddoutput(qf_diff,'dv',d,ceed_eval_grad,err) 169 170! Operator - apply 171 call ceedoperatorcreate(ceed,qf_diff,ceed_qfunction_none,& 172 & ceed_qfunction_none,op_diff,err) 173 call ceedoperatorsetfield(op_diff,'du',erestrictu,& 174 & bu,ceed_vector_active,err) 175 call ceedoperatorsetfield(op_diff,'qdata',erestrictqi,& 176 & ceed_basis_collocated,qdata,err) 177 call ceedoperatorsetfield(op_diff,'dv',erestrictu,& 178 & bu,ceed_vector_active,err) 179 180! Apply original Poisson Operator 181 call ceedvectorcreate(ceed,ndofs,u,err) 182 call ceedvectorsetvalue(u,1.d0,err) 183 call ceedvectorcreate(ceed,ndofs,v,err) 184 call ceedvectorsetvalue(v,0.d0,err) 185 call ceedoperatorapply(op_diff,u,v,ceed_request_immediate,err) 186 187! Check Output 188 call ceedvectorgetarrayread(v,ceed_mem_host,vv,voffset,err) 189 do i=1,ndofs 190 if (abs(vv(voffset+i))>1.0d-14) then 191! LCOV_EXCL_START 192 write(*,*) 'Error: Operator computed v[i] = ',vv(voffset+i),' != 0.0' 193! LCOV_EXCL_STOP 194 endif 195 enddo 196 call ceedvectorrestorearrayread(v,vv,voffset,err) 197 198! Assemble QFunction 199 call ceedoperatorassemblelinearqfunction(op_diff,a,erestrictlini,& 200 & ceed_request_immediate,err) 201 202! QFunction - apply linearized 203 call ceedqfunctioncreateinterior(ceed,1,diff_lin,& 204 &SOURCE_DIR& 205 &//'t531-operator.h:diff_lin'//char(0),qf_diff_lin,err) 206 call ceedqfunctionaddinput(qf_diff_lin,'du',d,ceed_eval_grad,err) 207 call ceedqfunctionaddinput(qf_diff_lin,'qdata',d*d,ceed_eval_none,err) 208 call ceedqfunctionaddoutput(qf_diff_lin,'dv',d,ceed_eval_grad,err) 209 210! Operator - apply linearized 211 call ceedoperatorcreate(ceed,qf_diff_lin,ceed_qfunction_none,& 212 & ceed_qfunction_none,op_diff_lin,err) 213 call ceedoperatorsetfield(op_diff_lin,'du',erestrictu,& 214 & bu,ceed_vector_active,err) 215 call ceedoperatorsetfield(op_diff_lin,'qdata',erestrictlini,& 216 & ceed_basis_collocated,a,err) 217 call ceedoperatorsetfield(op_diff_lin,'dv',erestrictu,& 218 & bu,ceed_vector_active,err) 219 220! Apply linearized Poisson Operator 221 call ceedvectorsetvalue(v,0.d0,err) 222 call ceedoperatorapply(op_diff_lin,u,v,ceed_request_immediate,err) 223 224! Check Output 225 call ceedvectorgetarrayread(v,ceed_mem_host,vv,voffset,err) 226 do i=1,ndofs 227 if (abs(vv(voffset+i))>1.0d-14) then 228! LCOV_EXCL_START 229 write(*,*) 'Error: Linearized operator computed v[i] = ',vv(voffset+i),& 230 & ' != 0.0' 231! LCOV_EXCL_STOP 232 endif 233 enddo 234 call ceedvectorrestorearrayread(v,vv,voffset,err) 235 236! Cleanup 237 call ceedqfunctiondestroy(qf_setup,err) 238 call ceedqfunctiondestroy(qf_diff,err) 239 call ceedqfunctiondestroy(qf_diff_lin,err) 240 call ceedoperatordestroy(op_setup,err) 241 call ceedoperatordestroy(op_diff,err) 242 call ceedoperatordestroy(op_diff_lin,err) 243 call ceedelemrestrictiondestroy(erestrictu,err) 244 call ceedelemrestrictiondestroy(erestrictx,err) 245 call ceedelemrestrictiondestroy(erestrictxi,err) 246 call ceedelemrestrictiondestroy(erestrictui,err) 247 call ceedelemrestrictiondestroy(erestrictqi,err) 248 call ceedelemrestrictiondestroy(erestrictlini,err) 249 call ceedbasisdestroy(bu,err) 250 call ceedbasisdestroy(bx,err) 251 call ceedvectordestroy(x,err) 252 call ceedvectordestroy(a,err) 253 call ceedvectordestroy(u,err) 254 call ceedvectordestroy(v,err) 255 call ceedvectordestroy(qdata,err) 256 call ceeddestroy(ceed,err) 257 end 258!----------------------------------------------------------------------- 259