Skip to content

Flux moments are not printed correctly when fluxp == 2 #23

@davidz-cerebras

Description

@davidz-cerebras

Setting fluxp == 2 is supposed to cause all flux moments to be output to the flux file. However, the Fortran reference implementation of SUBROUTINE output_flux_file erroneously prints multiple copies of moment 1 (the isotropic component) where moments 2 through cmom (the anisotropic components) are supposed to be.

I think the offending statement is here on line 332 of output.f90, which only reads from flux0 instead of conditionally reading from flux0 or fluxm depending on the value of l. Note that the output_send calls on line 322 do correctly implement this conditional logic:

            IF ( l == 1 ) THEN
              CALL output_send ( mtag, flux0(:,:,kloc,g) )
            ELSE
              CALL output_send ( mtag, fluxm(l-1,:,:,kloc,g) )
            END IF

However, in non-MPI runs, output_send is a no-op, causing the program to revert to faulty single-rank behavior. In addition, the FORMAT statement on line 370 uses the I1 specifier:

    323 FORMAT( /, 2X, 'Moment = ', I1 )

This causes the erroneous output Moment = * to be printed whenever the moment index is 10 or greater (which commonly occurs in 2D and 3D problems whenever nmom >= 4). The I0 specifier should be used instead to allow arbitrarily large integers to be printed without unnecessary padding.

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions