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