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