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