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 stridesu(3),stridesqd(3) 67 integer erestrictx,erestrictu,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 126 call ceedelemrestrictioncreate(ceed,imode,nelem,p*p,ndofs,1,& 127 & ceed_mem_host,ceed_use_pointer,indx,erestrictu,err) 128 stridesu=[1,q*q,q*q] 129 call ceedelemrestrictioncreatestrided(ceed,nelem,q*q,nqpts,& 130 & 1,stridesu,erestrictui,err) 131 132 stridesqd=[1,q*q,q*q*d*(d+1)/2] 133 call ceedelemrestrictioncreatestrided(ceed,nelem,q*q,nqpts,& 134 & d*(d+1)/2,stridesqd,erestrictqi,err) 135 136! Bases 137 call ceedbasiscreatetensorh1lagrange(ceed,d,d,p,q,ceed_gauss,& 138 & bx,err) 139 call ceedbasiscreatetensorh1lagrange(ceed,d,1,p,q,ceed_gauss,& 140 & bu,err) 141 142! QFunction - setup 143 call ceedqfunctioncreateinterior(ceed,1,setup,& 144 &SOURCE_DIR& 145 &//'t531-operator.h:setup'//char(0),qf_setup,err) 146 call ceedqfunctionaddinput(qf_setup,'dx',d*d,ceed_eval_grad,err) 147 call ceedqfunctionaddinput(qf_setup,'_weight',1,ceed_eval_weight,err) 148 call ceedqfunctionaddoutput(qf_setup,'qdata',d*(d+1)/2,ceed_eval_none,err) 149 150! Operator - setup 151 call ceedoperatorcreate(ceed,qf_setup,ceed_qfunction_none,& 152 & ceed_qfunction_none,op_setup,err) 153 call ceedoperatorsetfield(op_setup,'dx',erestrictx,& 154 & bx,ceed_vector_active,err) 155 call ceedoperatorsetfield(op_setup,'_weight',ceed_elemrestriction_none,& 156 & bx,ceed_vector_none,err) 157 call ceedoperatorsetfield(op_setup,'qdata',erestrictqi,& 158 & ceed_basis_collocated,ceed_vector_active,err) 159 160! Apply Setup Operator 161 call ceedoperatorapply(op_setup,x,qdata,ceed_request_immediate,err) 162 163! QFunction - apply 164 call ceedqfunctioncreateinterior(ceed,1,diff,& 165 &SOURCE_DIR& 166 &//'t531-operator.h:diff'//char(0),qf_diff,err) 167 call ceedqfunctionaddinput(qf_diff,'du',d,ceed_eval_grad,err) 168 call ceedqfunctionaddinput(qf_diff,'qdata',d*(d+1)/2,ceed_eval_none,err) 169 call ceedqfunctionaddoutput(qf_diff,'dv',d,ceed_eval_grad,err) 170 171! Operator - apply 172 call ceedoperatorcreate(ceed,qf_diff,ceed_qfunction_none,& 173 & ceed_qfunction_none,op_diff,err) 174 call ceedoperatorsetfield(op_diff,'du',erestrictu,& 175 & bu,ceed_vector_active,err) 176 call ceedoperatorsetfield(op_diff,'qdata',erestrictqi,& 177 & ceed_basis_collocated,qdata,err) 178 call ceedoperatorsetfield(op_diff,'dv',erestrictu,& 179 & bu,ceed_vector_active,err) 180 181! Apply original Poisson Operator 182 call ceedvectorcreate(ceed,ndofs,u,err) 183 call ceedvectorsetvalue(u,1.d0,err) 184 call ceedvectorcreate(ceed,ndofs,v,err) 185 call ceedvectorsetvalue(v,0.d0,err) 186 call ceedoperatorapply(op_diff,u,v,ceed_request_immediate,err) 187 188! Check Output 189 call ceedvectorgetarrayread(v,ceed_mem_host,vv,voffset,err) 190 do i=1,ndofs 191 if (abs(vv(voffset+i))>1.0d-14) then 192! LCOV_EXCL_START 193 write(*,*) 'Error: Operator computed v[i] = ',vv(voffset+i),' != 0.0' 194! LCOV_EXCL_STOP 195 endif 196 enddo 197 call ceedvectorrestorearrayread(v,vv,voffset,err) 198 199! Assemble QFunction 200 call ceedoperatorassemblelinearqfunction(op_diff,a,erestrictlini,& 201 & ceed_request_immediate,err) 202 203! QFunction - apply linearized 204 call ceedqfunctioncreateinterior(ceed,1,diff_lin,& 205 &SOURCE_DIR& 206 &//'t531-operator.h:diff_lin'//char(0),qf_diff_lin,err) 207 call ceedqfunctionaddinput(qf_diff_lin,'du',d,ceed_eval_grad,err) 208 call ceedqfunctionaddinput(qf_diff_lin,'qdata',d*d,ceed_eval_none,err) 209 call ceedqfunctionaddoutput(qf_diff_lin,'dv',d,ceed_eval_grad,err) 210 211! Operator - apply linearized 212 call ceedoperatorcreate(ceed,qf_diff_lin,ceed_qfunction_none,& 213 & ceed_qfunction_none,op_diff_lin,err) 214 call ceedoperatorsetfield(op_diff_lin,'du',erestrictu,& 215 & bu,ceed_vector_active,err) 216 call ceedoperatorsetfield(op_diff_lin,'qdata',erestrictlini,& 217 & ceed_basis_collocated,a,err) 218 call ceedoperatorsetfield(op_diff_lin,'dv',erestrictu,& 219 & bu,ceed_vector_active,err) 220 221! Apply linearized Poisson Operator 222 call ceedvectorsetvalue(v,0.d0,err) 223 call ceedoperatorapply(op_diff_lin,u,v,ceed_request_immediate,err) 224 225! Check Output 226 call ceedvectorgetarrayread(v,ceed_mem_host,vv,voffset,err) 227 do i=1,ndofs 228 if (abs(vv(voffset+i))>1.0d-14) then 229! LCOV_EXCL_START 230 write(*,*) 'Error: Linearized operator computed v[i] = ',vv(voffset+i),& 231 & ' != 0.0' 232! LCOV_EXCL_STOP 233 endif 234 enddo 235 call ceedvectorrestorearrayread(v,vv,voffset,err) 236 237! Cleanup 238 call ceedqfunctiondestroy(qf_setup,err) 239 call ceedqfunctiondestroy(qf_diff,err) 240 call ceedqfunctiondestroy(qf_diff_lin,err) 241 call ceedoperatordestroy(op_setup,err) 242 call ceedoperatordestroy(op_diff,err) 243 call ceedoperatordestroy(op_diff_lin,err) 244 call ceedelemrestrictiondestroy(erestrictu,err) 245 call ceedelemrestrictiondestroy(erestrictx,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