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(i)=u2(i)*u1(i) 27*d99fa3c5SJeremy L Thompson v1(q+i)=u2(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,bufine,bucoarse,bctof 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,ioffset 61*d99fa3c5SJeremy L Thompson real*8 val 62*d99fa3c5SJeremy L Thompson integer interpsize 63*d99fa3c5SJeremy L Thompson parameter(interpsize=pcoarse*pfine); 64*d99fa3c5SJeremy L Thompson real*8 interp1d(interpsize),interpctof(interpsize) 65*d99fa3c5SJeremy L Thompson 66*d99fa3c5SJeremy L Thompson real*8 hv(nufine*ncomp) 67*d99fa3c5SJeremy L Thompson real*8 total 68*d99fa3c5SJeremy L Thompson 69*d99fa3c5SJeremy L Thompson character arg*32 70*d99fa3c5SJeremy L Thompson 71*d99fa3c5SJeremy L Thompson external setup,mass 72*d99fa3c5SJeremy L Thompson 73*d99fa3c5SJeremy L Thompson call getarg(1,arg) 74*d99fa3c5SJeremy L Thompson call ceedinit(trim(arg)//char(0),ceed,err) 75*d99fa3c5SJeremy L Thompson 76*d99fa3c5SJeremy L Thompson do i=0,nx-1 77*d99fa3c5SJeremy L Thompson arrx(i+1)=i/(nx-1.d0) 78*d99fa3c5SJeremy L Thompson enddo 79*d99fa3c5SJeremy L Thompson do i=0,nelem-1 80*d99fa3c5SJeremy L Thompson indx(2*i+1)=i 81*d99fa3c5SJeremy L Thompson indx(2*i+2)=i+1 82*d99fa3c5SJeremy L Thompson enddo 83*d99fa3c5SJeremy L Thompson call ceedelemrestrictioncreate(ceed,nelem,2,1,1,nx,ceed_mem_host,& 84*d99fa3c5SJeremy L Thompson & ceed_use_pointer,indx,erestrictx,err) 85*d99fa3c5SJeremy L Thompson 86*d99fa3c5SJeremy L Thompson do i=0,nelem-1 87*d99fa3c5SJeremy L Thompson do j=0,pcoarse-1 88*d99fa3c5SJeremy L Thompson inducoarse(pcoarse*i+j+1)=i*(pcoarse-1)+j 89*d99fa3c5SJeremy L Thompson enddo 90*d99fa3c5SJeremy L Thompson enddo 91*d99fa3c5SJeremy L Thompson call ceedelemrestrictioncreate(ceed,nelem,pcoarse,ncomp,nucoarse,& 92*d99fa3c5SJeremy L Thompson & ncomp*nucoarse,ceed_mem_host,ceed_use_pointer,inducoarse,& 93*d99fa3c5SJeremy L Thompson & erestrictucoarse,err) 94*d99fa3c5SJeremy L Thompson 95*d99fa3c5SJeremy L Thompson do i=0,nelem-1 96*d99fa3c5SJeremy L Thompson do j=0,pfine-1 97*d99fa3c5SJeremy L Thompson indufine(pfine*i+j+1)=i*(pfine-1)+j 98*d99fa3c5SJeremy L Thompson enddo 99*d99fa3c5SJeremy L Thompson enddo 100*d99fa3c5SJeremy L Thompson call ceedelemrestrictioncreate(ceed,nelem,pfine,ncomp,nufine,& 101*d99fa3c5SJeremy L Thompson & ncomp*nufine,ceed_mem_host,ceed_use_pointer,indufine,& 102*d99fa3c5SJeremy L Thompson & erestrictufine,err) 103*d99fa3c5SJeremy L Thompson 104*d99fa3c5SJeremy L Thompson stridesu=[1,q,q] 105*d99fa3c5SJeremy L Thompson call ceedelemrestrictioncreatestrided(ceed,nelem,q,1,q*nelem,stridesu,& 106*d99fa3c5SJeremy L Thompson & erestrictui,err) 107*d99fa3c5SJeremy L Thompson 108*d99fa3c5SJeremy L Thompson call ceedbasiscreatetensorh1lagrange(ceed,1,1,2,q,ceed_gauss,bx,err) 109*d99fa3c5SJeremy L Thompson call ceedbasiscreatetensorh1lagrange(ceed,1,ncomp,pfine,q,ceed_gauss,& 110*d99fa3c5SJeremy L Thompson & bufine,err) 111*d99fa3c5SJeremy L Thompson 112*d99fa3c5SJeremy L Thompson call ceedqfunctioncreateinterior(ceed,1,setup,& 113*d99fa3c5SJeremy L Thompson &SOURCE_DIR& 114*d99fa3c5SJeremy L Thompson &//'t502-operator.h:setup'//char(0),qf_setup,err) 115*d99fa3c5SJeremy L Thompson call ceedqfunctionaddinput(qf_setup,'weights',1,ceed_eval_weight,err) 116*d99fa3c5SJeremy L Thompson call ceedqfunctionaddinput(qf_setup,'dx',1,ceed_eval_grad,err) 117*d99fa3c5SJeremy L Thompson call ceedqfunctionaddoutput(qf_setup,'qdata',1,ceed_eval_none,err) 118*d99fa3c5SJeremy L Thompson 119*d99fa3c5SJeremy L Thompson call ceedqfunctioncreateinterior(ceed,1,mass,& 120*d99fa3c5SJeremy L Thompson &SOURCE_DIR& 121*d99fa3c5SJeremy L Thompson &//'t502-operator.h:mass'//char(0),qf_mass,err) 122*d99fa3c5SJeremy L Thompson call ceedqfunctionaddinput(qf_mass,'qdata',1,ceed_eval_none,err) 123*d99fa3c5SJeremy L Thompson call ceedqfunctionaddinput(qf_mass,'u',ncomp,ceed_eval_interp,err) 124*d99fa3c5SJeremy L Thompson call ceedqfunctionaddoutput(qf_mass,'v',ncomp,ceed_eval_interp,err) 125*d99fa3c5SJeremy L Thompson 126*d99fa3c5SJeremy L Thompson call ceedoperatorcreate(ceed,qf_setup,ceed_qfunction_none,& 127*d99fa3c5SJeremy L Thompson & ceed_qfunction_none,op_setup,err) 128*d99fa3c5SJeremy L Thompson call ceedoperatorcreate(ceed,qf_mass,ceed_qfunction_none,& 129*d99fa3c5SJeremy L Thompson & ceed_qfunction_none,op_massfine,err) 130*d99fa3c5SJeremy L Thompson 131*d99fa3c5SJeremy L Thompson call ceedvectorcreate(ceed,nx,x,err) 132*d99fa3c5SJeremy L Thompson xoffset=0 133*d99fa3c5SJeremy L Thompson call ceedvectorsetarray(x,ceed_mem_host,ceed_use_pointer,arrx,xoffset,err) 134*d99fa3c5SJeremy L Thompson call ceedvectorcreate(ceed,nelem*q,qdata,err) 135*d99fa3c5SJeremy L Thompson 136*d99fa3c5SJeremy L Thompson call ceedoperatorsetfield(op_setup,'weights',ceed_elemrestriction_none,& 137*d99fa3c5SJeremy L Thompson & bx,ceed_vector_none,err) 138*d99fa3c5SJeremy L Thompson call ceedoperatorsetfield(op_setup,'dx',erestrictx,bx,& 139*d99fa3c5SJeremy L Thompson & ceed_vector_active,err) 140*d99fa3c5SJeremy L Thompson call ceedoperatorsetfield(op_setup,'qdata',erestrictui,& 141*d99fa3c5SJeremy L Thompson & ceed_basis_collocated,ceed_vector_active,err) 142*d99fa3c5SJeremy L Thompson call ceedoperatorsetfield(op_massfine,'qdata',erestrictui,& 143*d99fa3c5SJeremy L Thompson & ceed_basis_collocated,qdata,err) 144*d99fa3c5SJeremy L Thompson call ceedoperatorsetfield(op_massfine,'u',erestrictufine,bufine,& 145*d99fa3c5SJeremy L Thompson & ceed_vector_active,err) 146*d99fa3c5SJeremy L Thompson call ceedoperatorsetfield(op_massfine,'v',erestrictufine,bufine,& 147*d99fa3c5SJeremy L Thompson & ceed_vector_active,err) 148*d99fa3c5SJeremy L Thompson 149*d99fa3c5SJeremy L Thompson call ceedoperatorapply(op_setup,x,qdata,ceed_request_immediate,err) 150*d99fa3c5SJeremy L Thompson 151*d99fa3c5SJeremy L Thompson! Create multigrid level 152*d99fa3c5SJeremy L Thompson call ceedvectorcreate(ceed,ncomp*nufine,pmultfine,err) 153*d99fa3c5SJeremy L Thompson val=1.0 154*d99fa3c5SJeremy L Thompson call ceedvectorsetvalue(pmultfine,val,err) 155*d99fa3c5SJeremy L Thompson call ceedbasiscreatetensorh1lagrange(ceed,1,ncomp,pcoarse,q,& 156*d99fa3c5SJeremy L Thompson & ceed_gauss,bucoarse,err) 157*d99fa3c5SJeremy L Thompson call ceedbasiscreatetensorh1lagrange(ceed,1,ncomp,pcoarse,pfine,& 158*d99fa3c5SJeremy L Thompson & ceed_gauss_lobatto,bctof,err) 159*d99fa3c5SJeremy L Thompson call ceedbasisgetinterp1d(bctof,interp1d,ioffset,err) 160*d99fa3c5SJeremy L Thompson do i=1,interpsize 161*d99fa3c5SJeremy L Thompson interpctof(i)=interp1d(ioffset+i) 162*d99fa3c5SJeremy L Thompson enddo 163*d99fa3c5SJeremy L Thompson call ceedoperatormultigridlevelcreateh1(op_massfine,pmultfine,& 164*d99fa3c5SJeremy L Thompson & erestrictucoarse,bucoarse,interpctof,op_masscoarse,& 165*d99fa3c5SJeremy L Thompson & op_prolong,op_restrict,err) 166*d99fa3c5SJeremy L Thompson 167*d99fa3c5SJeremy L Thompson! Coarse problem 168*d99fa3c5SJeremy L Thompson call ceedvectorcreate(ceed,ncomp*nucoarse,ucoarse,err) 169*d99fa3c5SJeremy L Thompson val=1.0 170*d99fa3c5SJeremy L Thompson call ceedvectorsetvalue(ucoarse,val,err) 171*d99fa3c5SJeremy L Thompson call ceedvectorcreate(ceed,ncomp*nucoarse,vcoarse,err) 172*d99fa3c5SJeremy L Thompson call ceedoperatorapply(op_masscoarse,ucoarse,vcoarse,& 173*d99fa3c5SJeremy L Thompson & ceed_request_immediate,err) 174*d99fa3c5SJeremy L Thompson 175*d99fa3c5SJeremy L Thompson! Check output 176*d99fa3c5SJeremy L Thompson call ceedvectorgetarrayread(vcoarse,ceed_mem_host,hv,voffset,err) 177*d99fa3c5SJeremy L Thompson total=0. 178*d99fa3c5SJeremy L Thompson do i=1,nucoarse*ncomp 179*d99fa3c5SJeremy L Thompson total=total+hv(voffset+i) 180*d99fa3c5SJeremy L Thompson enddo 181*d99fa3c5SJeremy L Thompson if (abs(total-2.)>1.0d-10) then 182*d99fa3c5SJeremy L Thompson! LCOV_EXCL_START 183*d99fa3c5SJeremy L Thompson write(*,*) 'Computed Area: ',total,' != True Area: 1.0' 184*d99fa3c5SJeremy L Thompson! LCOV_EXCL_STOP 185*d99fa3c5SJeremy L Thompson endif 186*d99fa3c5SJeremy L Thompson call ceedvectorrestorearrayread(vcoarse,hv,voffset,err) 187*d99fa3c5SJeremy L Thompson 188*d99fa3c5SJeremy L Thompson! Prolong coarse u 189*d99fa3c5SJeremy L Thompson call ceedvectorcreate(ceed,ncomp*nufine,ufine,err) 190*d99fa3c5SJeremy L Thompson call ceedoperatorapply(op_prolong,ucoarse,ufine,& 191*d99fa3c5SJeremy L Thompson & ceed_request_immediate,err) 192*d99fa3c5SJeremy L Thompson 193*d99fa3c5SJeremy L Thompson! Fine problem 194*d99fa3c5SJeremy L Thompson call ceedvectorcreate(ceed,ncomp*nufine,vfine,err) 195*d99fa3c5SJeremy L Thompson call ceedoperatorapply(op_massfine,ufine,vfine,& 196*d99fa3c5SJeremy L Thompson & ceed_request_immediate,err) 197*d99fa3c5SJeremy L Thompson 198*d99fa3c5SJeremy L Thompson! Check output 199*d99fa3c5SJeremy L Thompson call ceedvectorgetarrayread(vfine,ceed_mem_host,hv,voffset,err) 200*d99fa3c5SJeremy L Thompson total=0. 201*d99fa3c5SJeremy L Thompson do i=1,nufine*ncomp 202*d99fa3c5SJeremy L Thompson total=total+hv(voffset+i) 203*d99fa3c5SJeremy L Thompson enddo 204*d99fa3c5SJeremy L Thompson if (abs(total-2.)>1.0d-10) then 205*d99fa3c5SJeremy L Thompson! LCOV_EXCL_START 206*d99fa3c5SJeremy L Thompson write(*,*) 'Computed Area Fine Grid: ',total,' != True Area: 1.0' 207*d99fa3c5SJeremy L Thompson! LCOV_EXCL_STOP 208*d99fa3c5SJeremy L Thompson endif 209*d99fa3c5SJeremy L Thompson call ceedvectorrestorearrayread(vfine,hv,voffset,err) 210*d99fa3c5SJeremy L Thompson 211*d99fa3c5SJeremy L Thompson! Restrict state to coarse grid 212*d99fa3c5SJeremy L Thompson call ceedoperatorapply(op_restrict,vfine,vcoarse,& 213*d99fa3c5SJeremy L Thompson & ceed_request_immediate,err) 214*d99fa3c5SJeremy L Thompson 215*d99fa3c5SJeremy L Thompson! Check output 216*d99fa3c5SJeremy L Thompson call ceedvectorgetarrayread(vcoarse,ceed_mem_host,hv,voffset,err) 217*d99fa3c5SJeremy L Thompson total=0. 218*d99fa3c5SJeremy L Thompson do i=1,nucoarse*ncomp 219*d99fa3c5SJeremy L Thompson total=total+hv(voffset+i) 220*d99fa3c5SJeremy L Thompson enddo 221*d99fa3c5SJeremy L Thompson if (abs(total-2.)>1.0d-10) then 222*d99fa3c5SJeremy L Thompson! LCOV_EXCL_START 223*d99fa3c5SJeremy L Thompson write(*,*) 'Computed Area Coarse Grid: ',total,' != True Area: 1.0' 224*d99fa3c5SJeremy L Thompson! LCOV_EXCL_STOP 225*d99fa3c5SJeremy L Thompson endif 226*d99fa3c5SJeremy L Thompson call ceedvectorrestorearrayread(vcoarse,hv,voffset,err) 227*d99fa3c5SJeremy L Thompson 228*d99fa3c5SJeremy L Thompson call ceedvectordestroy(qdata,err) 229*d99fa3c5SJeremy L Thompson call ceedvectordestroy(x,err) 230*d99fa3c5SJeremy L Thompson call ceedvectordestroy(ucoarse,err) 231*d99fa3c5SJeremy L Thompson call ceedvectordestroy(ufine,err) 232*d99fa3c5SJeremy L Thompson call ceedvectordestroy(vcoarse,err) 233*d99fa3c5SJeremy L Thompson call ceedvectordestroy(vfine,err) 234*d99fa3c5SJeremy L Thompson call ceedvectordestroy(pmultfine,err) 235*d99fa3c5SJeremy L Thompson call ceedoperatordestroy(op_masscoarse,err) 236*d99fa3c5SJeremy L Thompson call ceedoperatordestroy(op_massfine,err) 237*d99fa3c5SJeremy L Thompson call ceedoperatordestroy(op_prolong,err) 238*d99fa3c5SJeremy L Thompson call ceedoperatordestroy(op_restrict,err) 239*d99fa3c5SJeremy L Thompson call ceedoperatordestroy(op_setup,err) 240*d99fa3c5SJeremy L Thompson call ceedqfunctiondestroy(qf_mass,err) 241*d99fa3c5SJeremy L Thompson call ceedqfunctiondestroy(qf_setup,err) 242*d99fa3c5SJeremy L Thompson call ceedbasisdestroy(bufine,err) 243*d99fa3c5SJeremy L Thompson call ceedbasisdestroy(bx,err) 244*d99fa3c5SJeremy L Thompson call ceedelemrestrictiondestroy(erestrictucoarse,err) 245*d99fa3c5SJeremy L Thompson call ceedelemrestrictiondestroy(erestrictufine,err) 246*d99fa3c5SJeremy L Thompson call ceedelemrestrictiondestroy(erestrictx,err) 247*d99fa3c5SJeremy L Thompson call ceedelemrestrictiondestroy(erestrictui,err) 248*d99fa3c5SJeremy L Thompson call ceeddestroy(ceed,err) 249*d99fa3c5SJeremy L Thompson end 250*d99fa3c5SJeremy L Thompson!----------------------------------------------------------------------- 251