1!----------------------------------------------------------------------- 2! 3! Header with QFunctions 4! 5 include 't502-operator-f.h' 6!----------------------------------------------------------------------- 7 program test 8 9 include 'ceed/fortran.h' 10 11 integer ceed,err,i,j 12 integer stridesu(3) 13 integer erestrictx,erestrictui 14 integer erestrictucoarse,erestrictufine 15 integer bx,bufine,bucoarse,bctof 16 integer qf_setup,qf_mass 17 integer op_setup,op_masscoarse,op_massfine 18 integer op_prolong,op_restrict 19 integer qdata,x,ucoarse,ufine,vcoarse,vfine,pmultfine 20 integer nelem,pcoarse,pfine,q,ncomp 21 parameter(ncomp=2) 22 parameter(nelem=15) 23 parameter(pcoarse=3) 24 parameter(pfine=5) 25 parameter(q=8) 26 integer nx,nucoarse,nufine 27 parameter(nx=nelem+1) 28 parameter(nucoarse=nelem*(pcoarse-1)+1) 29 parameter(nufine=nelem*(pfine-1)+1) 30 integer indx(nelem*2) 31 integer inducoarse(nelem*pcoarse) 32 integer indufine(nelem*pfine) 33 real*8 arrx(nx) 34 integer*8 voffset,xoffset,ioffset 35 real*8 val 36 integer interpsize 37 parameter(interpsize=pcoarse*pfine); 38 real*8 interp1d(interpsize),interpctof(interpsize) 39 40 real*8 hv(nufine*ncomp) 41 real*8 total 42 43 character arg*32 44 45 external setup,mass 46 47 call getarg(1,arg) 48 call ceedinit(trim(arg)//char(0),ceed,err) 49 50 do i=0,nx-1 51 arrx(i+1)=i/(nx-1.d0) 52 enddo 53 do i=0,nelem-1 54 indx(2*i+1)=i 55 indx(2*i+2)=i+1 56 enddo 57 call ceedelemrestrictioncreate(ceed,nelem,2,1,1,nx,ceed_mem_host,& 58 & ceed_use_pointer,indx,erestrictx,err) 59 60 do i=0,nelem-1 61 do j=0,pcoarse-1 62 inducoarse(pcoarse*i+j+1)=i*(pcoarse-1)+j 63 enddo 64 enddo 65 call ceedelemrestrictioncreate(ceed,nelem,pcoarse,ncomp,nucoarse,& 66 & ncomp*nucoarse,ceed_mem_host,ceed_use_pointer,inducoarse,& 67 & erestrictucoarse,err) 68 69 do i=0,nelem-1 70 do j=0,pfine-1 71 indufine(pfine*i+j+1)=i*(pfine-1)+j 72 enddo 73 enddo 74 call ceedelemrestrictioncreate(ceed,nelem,pfine,ncomp,nufine,& 75 & ncomp*nufine,ceed_mem_host,ceed_use_pointer,indufine,& 76 & erestrictufine,err) 77 78 stridesu=[1,q,q] 79 call ceedelemrestrictioncreatestrided(ceed,nelem,q,1,q*nelem,stridesu,& 80 & erestrictui,err) 81 82 call ceedbasiscreatetensorh1lagrange(ceed,1,1,2,q,ceed_gauss,bx,err) 83 call ceedbasiscreatetensorh1lagrange(ceed,1,ncomp,pfine,q,ceed_gauss,& 84 & bufine,err) 85 86 call ceedqfunctioncreateinterior(ceed,1,setup,& 87 &SOURCE_DIR& 88 &//'t502-operator.h:setup'//char(0),qf_setup,err) 89 call ceedqfunctionaddinput(qf_setup,'weight',1,ceed_eval_weight,err) 90 call ceedqfunctionaddinput(qf_setup,'dx',1,ceed_eval_grad,err) 91 call ceedqfunctionaddoutput(qf_setup,'qdata',1,ceed_eval_none,err) 92 93 call ceedqfunctioncreateinterior(ceed,1,mass,& 94 &SOURCE_DIR& 95 &//'t502-operator.h:mass'//char(0),qf_mass,err) 96 call ceedqfunctionaddinput(qf_mass,'qdata',1,ceed_eval_none,err) 97 call ceedqfunctionaddinput(qf_mass,'u',ncomp,ceed_eval_interp,err) 98 call ceedqfunctionaddoutput(qf_mass,'v',ncomp,ceed_eval_interp,err) 99 100 call ceedoperatorcreate(ceed,qf_setup,ceed_qfunction_none,& 101 & ceed_qfunction_none,op_setup,err) 102 call ceedoperatorcreate(ceed,qf_mass,ceed_qfunction_none,& 103 & ceed_qfunction_none,op_massfine,err) 104 105 call ceedvectorcreate(ceed,nx,x,err) 106 xoffset=0 107 call ceedvectorsetarray(x,ceed_mem_host,ceed_use_pointer,arrx,xoffset,err) 108 call ceedvectorcreate(ceed,nelem*q,qdata,err) 109 110 call ceedoperatorsetfield(op_setup,'weight',ceed_elemrestriction_none,& 111 & bx,ceed_vector_none,err) 112 call ceedoperatorsetfield(op_setup,'dx',erestrictx,bx,& 113 & ceed_vector_active,err) 114 call ceedoperatorsetfield(op_setup,'qdata',erestrictui,& 115 ceed_basis_none,ceed_vector_active,err) 116 call ceedoperatorsetfield(op_massfine,'qdata',erestrictui,& 117 ceed_basis_none,qdata,err) 118 call ceedoperatorsetfield(op_massfine,'u',erestrictufine,bufine,& 119 & ceed_vector_active,err) 120 call ceedoperatorsetfield(op_massfine,'v',erestrictufine,bufine,& 121 & ceed_vector_active,err) 122 123 call ceedoperatorapply(op_setup,x,qdata,ceed_request_immediate,err) 124 125! Create multigrid level 126 call ceedvectorcreate(ceed,ncomp*nufine,pmultfine,err) 127 val=1.0 128 call ceedvectorsetvalue(pmultfine,val,err) 129 call ceedbasiscreatetensorh1lagrange(ceed,1,ncomp,pcoarse,q,& 130 & ceed_gauss,bucoarse,err) 131 call ceedbasiscreatetensorh1lagrange(ceed,1,ncomp,pcoarse,pfine,& 132 & ceed_gauss_lobatto,bctof,err) 133 call ceedbasisgetinterp1d(bctof,interp1d,ioffset,err) 134 do i=1,interpsize 135 interpctof(i)=interp1d(ioffset+i) 136 enddo 137 call ceedoperatormultigridlevelcreateh1(op_massfine,pmultfine,& 138 & erestrictucoarse,bucoarse,interpctof,op_masscoarse,& 139 & op_prolong,op_restrict,err) 140 141! Coarse problem 142 call ceedvectorcreate(ceed,ncomp*nucoarse,ucoarse,err) 143 val=1.0 144 call ceedvectorsetvalue(ucoarse,val,err) 145 call ceedvectorcreate(ceed,ncomp*nucoarse,vcoarse,err) 146 call ceedoperatorapply(op_masscoarse,ucoarse,vcoarse,& 147 & ceed_request_immediate,err) 148 149! Check output 150 call ceedvectorgetarrayread(vcoarse,ceed_mem_host,hv,voffset,err) 151 total=0. 152 do i=1,nucoarse*ncomp 153 total=total+hv(voffset+i) 154 enddo 155 if (abs(total-2.)>1.0d-10) then 156! LCOV_EXCL_START 157 write(*,*) 'Computed Area: ',total,' != True Area: 1.0' 158! LCOV_EXCL_STOP 159 endif 160 call ceedvectorrestorearrayread(vcoarse,hv,voffset,err) 161 162! Prolong coarse u 163 call ceedvectorcreate(ceed,ncomp*nufine,ufine,err) 164 call ceedoperatorapply(op_prolong,ucoarse,ufine,& 165 & ceed_request_immediate,err) 166 167! Fine problem 168 call ceedvectorcreate(ceed,ncomp*nufine,vfine,err) 169 call ceedoperatorapply(op_massfine,ufine,vfine,& 170 & ceed_request_immediate,err) 171 172! Check output 173 call ceedvectorgetarrayread(vfine,ceed_mem_host,hv,voffset,err) 174 total=0. 175 do i=1,nufine*ncomp 176 total=total+hv(voffset+i) 177 enddo 178 if (abs(total-2.)>1.0d-10) then 179! LCOV_EXCL_START 180 write(*,*) 'Computed Area Fine Grid: ',total,' != True Area: 1.0' 181! LCOV_EXCL_STOP 182 endif 183 call ceedvectorrestorearrayread(vfine,hv,voffset,err) 184 185! Restrict state to coarse grid 186 call ceedoperatorapply(op_restrict,vfine,vcoarse,& 187 & ceed_request_immediate,err) 188 189! Check output 190 call ceedvectorgetarrayread(vcoarse,ceed_mem_host,hv,voffset,err) 191 total=0. 192 do i=1,nucoarse*ncomp 193 total=total+hv(voffset+i) 194 enddo 195 if (abs(total-2.)>1.0d-10) then 196! LCOV_EXCL_START 197 write(*,*) 'Computed Area Coarse Grid: ',total,' != True Area: 1.0' 198! LCOV_EXCL_STOP 199 endif 200 call ceedvectorrestorearrayread(vcoarse,hv,voffset,err) 201 202 call ceedvectordestroy(qdata,err) 203 call ceedvectordestroy(x,err) 204 call ceedvectordestroy(ucoarse,err) 205 call ceedvectordestroy(ufine,err) 206 call ceedvectordestroy(vcoarse,err) 207 call ceedvectordestroy(vfine,err) 208 call ceedvectordestroy(pmultfine,err) 209 call ceedoperatordestroy(op_masscoarse,err) 210 call ceedoperatordestroy(op_massfine,err) 211 call ceedoperatordestroy(op_prolong,err) 212 call ceedoperatordestroy(op_restrict,err) 213 call ceedoperatordestroy(op_setup,err) 214 call ceedqfunctiondestroy(qf_mass,err) 215 call ceedqfunctiondestroy(qf_setup,err) 216 call ceedbasisdestroy(bufine,err) 217 call ceedbasisdestroy(bx,err) 218 call ceedelemrestrictiondestroy(erestrictucoarse,err) 219 call ceedelemrestrictiondestroy(erestrictufine,err) 220 call ceedelemrestrictiondestroy(erestrictx,err) 221 call ceedelemrestrictiondestroy(erestrictui,err) 222 call ceeddestroy(ceed,err) 223 end 224!----------------------------------------------------------------------- 225