xref: /libCEED/tests/t502-operator-f.f90 (revision 39b2de37682296be8460181179eb4e44de5cc3de)
1!-----------------------------------------------------------------------
2      subroutine setup(ctx,q,u1,u2,u3,u4,u5,u6,u7,u8,u9,u10,u11,u12,u13,u14,&
3&           u15,u16,v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,v13,v14,v15,v16,ierr)
4      real*8 ctx
5      real*8 u1(1)
6      real*8 u2(1)
7      real*8 v1(1)
8      integer q,ierr
9
10      do i=1,q
11        v1(i)=u1(i)*u2(i)
12      enddo
13
14      ierr=0
15      end
16!-----------------------------------------------------------------------
17      subroutine mass(ctx,q,u1,u2,u3,u4,u5,u6,u7,u8,u9,u10,u11,u12,u13,u14,&
18&           u15,u16,v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,v13,v14,v15,v16,ierr)
19      real*8 ctx
20      real*8 u1(1)
21      real*8 u2(1)
22      real*8 v1(1)
23      integer q,ierr
24
25      do i=1,q
26        v1(i)=u2(i)*u1(i)
27        v1(q+i)=u2(q+i)*u1(i)
28      enddo
29
30      ierr=0
31      end
32!-----------------------------------------------------------------------
33      program test
34
35      include 'ceedf.h'
36
37      integer ceed,err,i,j
38      integer erestrictx,erestrictu,erestrictxi,erestrictui
39      integer bx,bu
40      integer qf_setup,qf_mass
41      integer op_setup,op_mass
42      integer qdata,x,u,v
43      integer nelem,p,q
44      parameter(nelem=15)
45      parameter(p=5)
46      parameter(q=8)
47      integer nx,nu
48      parameter(nx=nelem+1)
49      parameter(nu=nelem*(p-1)+1)
50      integer indx(nelem*2)
51      integer indu(nelem*p)
52      real*8 arrx(nx)
53      integer*8 voffset,xoffset
54
55      real*8 hu(nu*2),hv(nu*2)
56      real*8 total1,total2
57
58      character arg*32
59
60      external setup,mass
61
62      call getarg(1,arg)
63      call ceedinit(trim(arg)//char(0),ceed,err)
64
65      do i=0,nx-1
66        arrx(i+1)=i/(nx-1.d0)
67      enddo
68      do i=0,nelem-1
69        indx(2*i+1)=i
70        indx(2*i+2)=i+1
71      enddo
72
73      call ceedelemrestrictioncreate(ceed,nelem,2,nx,1,ceed_mem_host,&
74     & ceed_use_pointer,indx,erestrictx,err)
75      call ceedelemrestrictioncreateidentity(ceed,nelem,2,2*nelem,1,&
76     & erestrictxi,err)
77
78      do i=0,nelem-1
79        do j=0,p-1
80          indu(p*i+j+1)=i*(p-1)+j
81        enddo
82      enddo
83
84      call ceedelemrestrictioncreate(ceed,nelem,p,nu,2,ceed_mem_host,&
85     & ceed_use_pointer,indu,erestrictu,err)
86      call ceedelemrestrictioncreateidentity(ceed,nelem,q,q*nelem,1,&
87     & erestrictui,err)
88
89      call ceedbasiscreatetensorh1lagrange(ceed,1,1,2,q,ceed_gauss,bx,err)
90      call ceedbasiscreatetensorh1lagrange(ceed,1,2,p,q,ceed_gauss,bu,err)
91
92      call ceedqfunctioncreateinterior(ceed,1,setup,&
93     &SOURCE_DIR&
94     &//'t502-operator.h:setup'//char(0),qf_setup,err)
95      call ceedqfunctionaddinput(qf_setup,'_weight',1,ceed_eval_weight,err)
96      call ceedqfunctionaddinput(qf_setup,'dx',1,ceed_eval_grad,err)
97      call ceedqfunctionaddoutput(qf_setup,'rho',1,ceed_eval_none,err)
98
99      call ceedqfunctioncreateinterior(ceed,1,mass,&
100     &SOURCE_DIR&
101     &//'t502-operator.h:mass'//char(0),qf_mass,err)
102      call ceedqfunctionaddinput(qf_mass,'rho',1,ceed_eval_none,err)
103      call ceedqfunctionaddinput(qf_mass,'u',2,ceed_eval_interp,err)
104      call ceedqfunctionaddoutput(qf_mass,'v',2,ceed_eval_interp,err)
105
106      call ceedoperatorcreate(ceed,qf_setup,ceed_qfunction_none,&
107     & ceed_qfunction_none,op_setup,err)
108      call ceedoperatorcreate(ceed,qf_mass,ceed_qfunction_none,&
109     & ceed_qfunction_none,op_mass,err)
110
111      call ceedvectorcreate(ceed,nx,x,err)
112      xoffset=0
113      call ceedvectorsetarray(x,ceed_mem_host,ceed_use_pointer,arrx,xoffset,err)
114      call ceedvectorcreate(ceed,nelem*q,qdata,err)
115
116      call ceedoperatorsetfield(op_setup,'_weight',erestrictxi,&
117     & ceed_notranspose,bx,ceed_vector_none,err)
118      call ceedoperatorsetfield(op_setup,'dx',erestrictx,&
119     & ceed_notranspose,bx,ceed_vector_active,err)
120      call ceedoperatorsetfield(op_setup,'rho',erestrictui,&
121     & ceed_notranspose,ceed_basis_collocated,ceed_vector_active,err)
122      call ceedoperatorsetfield(op_mass,'rho',erestrictui,&
123     & ceed_notranspose,ceed_basis_collocated,qdata,err)
124      call ceedoperatorsetfield(op_mass,'u',erestrictu,&
125     & ceed_transpose,bu,ceed_vector_active,err)
126      call ceedoperatorsetfield(op_mass,'v',erestrictu,&
127     & ceed_transpose,bu,ceed_vector_active,err)
128
129      call ceedoperatorapply(op_setup,x,qdata,ceed_request_immediate,err)
130
131      call ceedvectorcreate(ceed,2*nu,u,err)
132      call ceedvectorgetarray(u,ceed_mem_host,hu,voffset,err)
133      do i=1,nu
134        hu(voffset+2*i-1)=1.
135        hu(voffset+2*i)=2.
136      enddo
137      call ceedvectorrestorearray(u,hu,voffset,err)
138      call ceedvectorcreate(ceed,2*nu,v,err)
139      call ceedoperatorapply(op_mass,u,v,ceed_request_immediate,err)
140
141      call ceedvectorgetarrayread(v,ceed_mem_host,hv,voffset,err)
142      total1=0.
143      total2=0.
144      do i=1,nu
145        total1=total1+hv(voffset+2*i-1)
146        total2=total2+hv(voffset+2*i)
147      enddo
148      if (abs(total1-1.)>1.0d-10) then
149! LCOV_EXCL_START
150        write(*,*) 'Computed Area: ',total1,' != True Area: 1.0'
151! LCOV_EXCL_STOP
152      endif
153      if (abs(total2-2.)>1.0d-10) then
154! LCOV_EXCL_START
155        write(*,*) 'Computed Area: ',total2,' != True Area: 2.0'
156! LCOV_EXCL_STOP
157      endif
158      call ceedvectorrestorearrayread(v,hv,voffset,err)
159
160      call ceedvectordestroy(qdata,err)
161      call ceedvectordestroy(x,err)
162      call ceedvectordestroy(u,err)
163      call ceedvectordestroy(v,err)
164      call ceedoperatordestroy(op_mass,err)
165      call ceedoperatordestroy(op_setup,err)
166      call ceedqfunctiondestroy(qf_mass,err)
167      call ceedqfunctiondestroy(qf_setup,err)
168      call ceedbasisdestroy(bu,err)
169      call ceedbasisdestroy(bx,err)
170      call ceedelemrestrictiondestroy(erestrictu,err)
171      call ceedelemrestrictiondestroy(erestrictx,err)
172      call ceedelemrestrictiondestroy(erestrictui,err)
173      call ceedelemrestrictiondestroy(erestrictxi,err)
174      call ceeddestroy(ceed,err)
175      end
176!-----------------------------------------------------------------------
177