1!----------------------------------------------------------------------- 2! 3! Header with QFunctions 4! 5 include 't510-operator-f.h' 6!----------------------------------------------------------------------- 7 program test 8 implicit none 9 include 'ceed/fortran.h' 10 11 integer ceed,err,i,j,k 12 integer stridesu(3) 13 integer erestrictx,erestrictu,erestrictui 14 integer bx,bu 15 integer qf_setup,qf_mass 16 integer op_setup,op_mass 17 integer qdata,x,a,u,v 18 integer nelem,p,q,d 19 integer row,col,offset 20 parameter(nelem=6) 21 parameter(p=3) 22 parameter(q=4) 23 parameter(d=2) 24 integer ndofs,nqpts,nx,ny 25 parameter(nx=3) 26 parameter(ny=2) 27 parameter(ndofs=(nx*2+1)*(ny*2+1)) 28 parameter(nqpts=nelem*q*q) 29 integer indx(nelem*p*p) 30 real*8 arrx(d*ndofs),aa(nqpts),uu(ndofs),vv(ndofs),atrue(ndofs) 31 integer*8 xoffset,aoffset,uoffset,voffset 32 33 character arg*32 34 35 external setup,mass 36 37 call getarg(1,arg) 38 39 call ceedinit(trim(arg)//char(0),ceed,err) 40 41! DoF Coordinates 42 do i=0,nx*2 43 do j=0,ny*2 44 arrx(i+j*(nx*2+1)+0*ndofs+1)=1.d0*i/(2*nx) 45 arrx(i+j*(nx*2+1)+1*ndofs+1)=1.d0*j/(2*ny) 46 enddo 47 enddo 48 call ceedvectorcreate(ceed,d*ndofs,x,err) 49 xoffset=0 50 call ceedvectorsetarray(x,ceed_mem_host,ceed_use_pointer,arrx,xoffset,err) 51 52! Qdata Vector 53 call ceedvectorcreate(ceed,nqpts,qdata,err) 54 55! Element Setup 56 do i=0,nelem-1 57 col=mod(i,nx) 58 row=i/nx 59 offset=col*(p-1)+row*(nx*2+1)*(p-1) 60 do j=0,p-1 61 do k=0,p-1 62 indx(p*(p*i+k)+j+1)=offset+k*(nx*2+1)+j 63 enddo 64 enddo 65 enddo 66 67! Restrictions 68 call ceedelemrestrictioncreate(ceed,nelem,p*p,d,ndofs,d*ndofs,& 69 & ceed_mem_host,ceed_use_pointer,indx,erestrictx,err) 70 71 call ceedelemrestrictioncreate(ceed,nelem,p*p,1,1,ndofs,& 72 & ceed_mem_host,ceed_use_pointer,indx,erestrictu,err) 73 stridesu=[1,q*q,q*q] 74 call ceedelemrestrictioncreatestrided(ceed,nelem,q*q,1,nqpts,& 75 & stridesu,erestrictui,err) 76 77! Bases 78 call ceedbasiscreatetensorh1lagrange(ceed,d,d,p,q,ceed_gauss,& 79 & bx,err) 80 call ceedbasiscreatetensorh1lagrange(ceed,d,1,p,q,ceed_gauss,& 81 & bu,err) 82 83! QFunctions 84! -- Setup 85 call ceedqfunctioncreateinterior(ceed,1,setup,& 86 &SOURCE_DIR& 87 &//'t510-operator.h:setup'//char(0),qf_setup,err) 88 call ceedqfunctionaddinput(qf_setup,'_weight',1,ceed_eval_weight,err) 89 call ceedqfunctionaddinput(qf_setup,'dx',d*d,ceed_eval_grad,err) 90 call ceedqfunctionaddoutput(qf_setup,'rho',1,ceed_eval_none,err) 91! -- Mass 92 call ceedqfunctioncreateinterior(ceed,1,mass,& 93 &SOURCE_DIR& 94 &//'t510-operator.h:mass'//char(0),qf_mass,err) 95 call ceedqfunctionaddinput(qf_mass,'rho',1,ceed_eval_none,err) 96 call ceedqfunctionaddinput(qf_mass,'u',1,ceed_eval_interp,err) 97 call ceedqfunctionaddoutput(qf_mass,'v',1,ceed_eval_interp,err) 98 99! Operators 100! -- Setup 101 call ceedoperatorcreate(ceed,qf_setup,ceed_qfunction_none,& 102 & ceed_qfunction_none,op_setup,err) 103 call ceedoperatorsetfield(op_setup,'_weight',ceed_elemrestriction_none,& 104 & bx,ceed_vector_none,err) 105 call ceedoperatorsetfield(op_setup,'dx',erestrictx,& 106 & bx,ceed_vector_active,err) 107 call ceedoperatorsetfield(op_setup,'rho',erestrictui,& 108 & ceed_basis_collocated,ceed_vector_active,err) 109! -- Mass 110 call ceedoperatorcreate(ceed,qf_mass,ceed_qfunction_none,& 111 & ceed_qfunction_none,op_mass,err) 112 call ceedoperatorsetfield(op_mass,'rho',erestrictui,& 113 & ceed_basis_collocated,qdata,err) 114 call ceedoperatorsetfield(op_mass,'u',erestrictu,& 115 & bu,ceed_vector_active,err) 116 call ceedoperatorsetfield(op_mass,'v',erestrictu,& 117 & bu,ceed_vector_active,err) 118 119! Apply Setup Operator 120 call ceedoperatorapply(op_setup,x,qdata,ceed_request_immediate,err) 121 122! Assemble Diagonal 123 call ceedvectorcreate(ceed,ndofs,a,err) 124 call ceedoperatorlinearassemblediagonal(op_mass,a,& 125 & ceed_request_immediate,err) 126 127! Manually assemble diagonal 128 call ceedvectorcreate(ceed,ndofs,u,err) 129 call ceedvectorsetvalue(u,0.d0,err) 130 call ceedvectorcreate(ceed,ndofs,v,err) 131 do i=1,ndofs 132 call ceedvectorgetarray(u,ceed_mem_host,uu,uoffset,err) 133 uu(i+uoffset)=1.d0 134 if (i>1) then 135 uu(i-1+uoffset)=0.d0 136 endif 137 call ceedvectorrestorearray(u,uu,uoffset,err) 138 139 call ceedoperatorapply(op_mass,u,v,ceed_request_immediate,err) 140 141 call ceedvectorgetarrayread(v,ceed_mem_host,vv,voffset,err) 142 atrue(i)=vv(voffset+i) 143 call ceedvectorrestorearrayread(v,vv,voffset,err) 144 enddo 145 146! Check Output 147 call ceedvectorgetarrayread(a,ceed_mem_host,aa,aoffset,err) 148 do i=1,ndofs 149 if (abs(aa(aoffset+i)-atrue(i))>1.0d-14) then 150! LCOV_EXCL_START 151 write(*,*) '[',i,'] Error in assembly: ',aa(aoffset+i),' != ',& 152 & atrue(i) 153! LCOV_EXCL_STOP 154 endif 155 enddo 156 call ceedvectorrestorearrayread(a,aa,aoffset,err) 157 158! Cleanup 159 call ceedqfunctiondestroy(qf_setup,err) 160 call ceedqfunctiondestroy(qf_mass,err) 161 call ceedoperatordestroy(op_setup,err) 162 call ceedoperatordestroy(op_mass,err) 163 call ceedelemrestrictiondestroy(erestrictu,err) 164 call ceedelemrestrictiondestroy(erestrictx,err) 165 call ceedelemrestrictiondestroy(erestrictui,err) 166 call ceedbasisdestroy(bu,err) 167 call ceedbasisdestroy(bx,err) 168 call ceedvectordestroy(x,err) 169 call ceedvectordestroy(a,err) 170 call ceedvectordestroy(u,err) 171 call ceedvectordestroy(v,err) 172 call ceedvectordestroy(qdata,err) 173 call ceeddestroy(ceed,err) 174 end 175!----------------------------------------------------------------------- 176