[LFRic] [PSyAD] PSyAD no longer computes adjoint of vector assigments properly in PSyclone 3.0-dev

[DrTVockerodtMO]

In testing PSyclone 3.0, we find that forward code that looks like

vector(:) = vector(:) + constant(:) * scalar

where vector and scalar are active no longer adjoints to

do idx = vector_ubound, vector_lbound, -1
  scalar = scalar + constant(idx) * vector(idx)
end do

but instead does

scalar = scalar + constant(:) * vector(:)

which fails to compile. The former behaviour is expected.

[arporter]

Thanks @DrTVockerodtMO, could you provide some Fortran source that shows this problem when running PSyAD?

[DrTVockerodtMO]

Thanks @DrTVockerodtMO, could you provide some Fortran source that shows this problem when running PSyAD?

Sure! For the forward strong_curl_kernel we have the lines

do df1 = 1,ndf1
  curl_u(:) = curl_u(:) + u(map1(df1)+k)*diff_basis_w1(:,df1,df2)
end do

these used to get adjointed pre 3.0 to

do df1 = ndf1, 1, -1
 do i = 3, 1, -1
   u(k + map1(df1)) = u(k + map1(df1)) + diff_basis_w1(i,df1,df2) * curl_u(i)
 enddo
enddo

which is correct, but now they get adjointed to

do df1 = ndf1, 1, -1
  u(k + map1(df1)) = u(k + map1(df1)) + diff_basis_w1(:,df1,df2) * curl_u(:)
enddo

which is not.

The other kernel that exhibit this new behaviour is vorticity_advection_kernel.

[arporter]

In the psyad examples, we do:

psyad -api lfric tl_vorticity_advection_kernel_mod_tweaked.F90 -oad ../adjoint_partial/adj_vorticity_advection_kernel_mod.F90 -a r_u res_dot_product vorticity_term cross_product1 cross_product2 j_vorticity u_at_quad mul2 vorticity_at_quad wind vorticity

and the problem is presumably to do with the ArrayAssignment2LoopsTrans. In fact, there seems to be two problems: that transformation appears to be failing and psyad is carrying on regardless of that failure.