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