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