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