stfc / PSyclone

Domain-specific compiler and code transformation system for Finite Difference/Volume/Element Earth-system models in Fortran
BSD 3-Clause "New" or "Revised" License
107 stars 28 forks source link

[psyad] Problem with adjoint of assignments involving scalars and array ranges #1581

Closed arporter closed 2 years ago

arporter commented 2 years ago
    do df = 1, ndf_wt
      theta_v_at_quad   = theta_v_at_quad                                 &
                        + theta_v_e(df)*wt_basis(1,df,qp1,qp2)
      grad_theta_v_at_quad(:) = grad_theta_v_at_quad(:)                   &
                               + theta_v_e(df)*wt_diff_basis(:,df,qp1,qp2)
    end do

gets transformed to:

      do df = ndf_wt, 1, -1
        theta_v_e(df) = theta_v_e(df) + grad_theta_v_at_quad(:) * wt_diff_basis(:,df,qp1,qp2)
        theta_v_e(df) = theta_v_e(df) + theta_v_at_quad * wt_basis(1,df,qp1,qp2)
      enddo

where the first assignment is not valid.

arporter commented 2 years ago

Really we have something like:

do ipt = 1, nqp_pts
  grad_theta_v_at_quad(ipt) = grad_theta_v_at_quad(ipt) + theta_v_e(df)*wt_diff_basis(ipt, df, qp1, qp2)
end do

grad_theta_v_at_quad and theta_v_e are both active. This should then become:

      do idx = 3, 1, -1
          theta_v_e(df) = theta_v_e(df) + grad_theta_v_at_quad(idx) * wt_diff_basis(idx,df,qp1,qp2)
      enddo

so we actually end up accumulating into theta_v_e(df).

Perhaps the easiest solution is to convert all implicit loops into explicit ones first @rupertford? (We can build on the NEMO transformation that does this.)

rupertford commented 2 years ago

Yes, I agree @arporter. We should be able to simply add the arrayrange2loop transformation to our internal pre-processing script.

The invalid Fortran that is produced is actually correct if you assume that subexpressions on the rhs can be implicit loops :-)

rupertford commented 2 years ago

Created branch 1581_psyad_array_range_to_loop