Getting a weird runtime error when trying to use Coarrays

Help
1

Hi all,
I am trying to use coarrays on a small piece of code I wrote earlier
Here is the coarray code (so far)

PROGRAM bmsolve_optimized_coarrays
  IMPLICIT NONE

  !Compiler options: -march=skylake -mtune=skylake -O3 -mfpmath=sse

  !Data type to store details of the beam cross-section
  TYPE data
     REAL, ALLOCATABLE :: disp(:), moment(:)
     REAL :: maximum_disp, maximum_values_disp(2)
     REAL :: pa=1.
     REAL :: pb,pc
     REAL :: In, A, d1, d2, d3, d4
     REAL :: fa, fb, fc
     REAL :: maximum_normal_stress, maximum_moment, maximum_values_moment(2)
     INTEGER :: id
  END type data

  !Parameters
  REAL, PARAMETER :: q=8.,disp_all=5.,pi=4.0*ATAN(1.0),l=6000,E=2E5
  INTEGER, PARAMETER :: lsize=6001, dmax=500

  !variables
  INTEGER :: i,j,k,r, maximum_disp_loc, maximum_disp_loc_coarray
  INTEGER(KIND=8) :: total_carray[*]
  REAL, ALLOCATABLE :: div(:),divl(:)
  REAL :: t1,t2
  TYPE(data) :: s
  TYPE(data) :: best_minimum_disp_carray[*],best_minimum_disp,best_arr(5)[*],best


  !Allocating arrays
  ALLOCATE(div(dmax));ALLOCATE(divl(lsize));ALLOCATE(s%disp(lsize));ALLOCATE(s%moment(lsize))


  !Output format
  100 FORMAT (X,A,25X,A)
  110 FORMAT (X,A,G23.8,3X,A)
  120 FORMAT (X,A,G36.8,3X,A)
  130 FORMAT (X,A,G25.8,2X,A)
  140 FORMAT (X,A,G27.8,X,A)
  150 FORMAT (X,A,G26.8,2X,A)
  160 FORMAT (X,A,G30.8,X,A)
  170 FORMAT (X,A,G26.8,2X,A)
  180 FORMAT (X,A,G18.8,1X,A)
  190 FORMAT (X,A,G20.8,X,A)
  191 FORMAT (X,A,G20.8,X,A)
  192 FORMAT (X,A,G33.8,X,A)

  DO i=1,lsize
     divl(i)=REAL(i)
  END DO

  DO i=1,dmax
     div(i)=REAL(i)
  END DO

  s%pb=2000.
  s%pc=4000.

  CALL find_react(s%fa, s%fb, s%fc, s%pb, s%pc)

  SYNC ALL()

  IF (THIS_IMAGE() .EQ. 1) THEN
     PRINT *, " "
     PRINT 192, "Initial reaction at A: ",s%fa,"N"
     PRINT 192, "Initial reaction at B: ",s%fb,"N"
     PRINT 192, "Initial reaction at C: ",s%fc,"N"
     PRINT *, " "
  END IF

  s%id=1
  s%d1=div(80)

  CALL find_area_I(s%id, s%d1, s%d2, s%d3, s%d4, s%A, S%In)
  CALL find_disp_moment(s%pb, s%pc, s%fb, s%fc, s%In, s%disp, s%moment)
  CALL find_maximum_disp_moment(s%disp, s%moment, s%maximum_values_disp, s%maximum_values_moment, s%maximum_disp, s%maximum_moment)

  IF (MAXLOC(s%maximum_values_disp,1)==1) THEN
     maximum_disp_loc=MAXLOC(s%disp,1)
  ELSE IF (MAXLOC(s%maximum_values_disp,1)==2) THEN
     maximum_disp_loc=MINLOC(s%disp,1)
  END IF

  IF (THIS_IMAGE() .EQ. 1) THEN
    PRINT 192, "Location of maximum displacement: ",(maximum_disp_loc-1),"mm"
    PRINT *, " "
  END IF!

  !Solid circular cross-section
  s%id=1
  best_arr(1)%maximum_disp=0.
  s%disp=0.
  s%d2=0.
  s%d3=0.
  s%d4=0.
  
  SYNC ALL()




 DEALLOCATE(div); DEALLOCATE(divl); DEALLOCATE(s%disp); DEALLOCATE(s%moment)

CONTAINS

  SUBROUTINE find_react(fa,fb,fc,pb,pc)
    IMPLICIT NONE

    !Data from main program
    REAL,INTENT(IN) :: pb,pc
    REAL,INTENT(OUT) :: fa,fb,fc
                                                      
    !Local variables
    REAL :: detA
    REAL :: A(2,2),A_in(2,2),X(2,1),B(2,1)


    A(1,1)=((pb**2)*((3*pb)-pb))/6.
    A(1,2)=((pb**2)*((3*pc)-pb))/6.
    A(2,1)=((pb**2)*((3*pc)-pb))/6.
    A(2,2)=((pc**2)*((3*pc)-pc))/6.

    B(1,1)=(q*(pb**2)*((6*(l**2))-(4*l*pb)+(pb**2)))/24.
    B(2,1)=(q*(pc**2)*((6*(l**2))-(4*l*pc)+(pc**2)))/24.

    detA=1/((A(1,1)*A(2,2))-(A(1,2)*A(2,1)))

    A_in(1,1)=+detA*A(2,2)
    A_in(2,1)=-detA*A(2,1)
    A_in(1,2)=-detA*A(1,2)
    A_in(2,2)=+detA*A(1,1)

    X(1,1)=(B(1,1)*A_in(1,1))+(A_in(1,2)*B(2,1))
    X(2,1)=(B(1,1)*A_in(2,1))+(A_in(2,2)*B(2,1))


    fb=X(1,1)
    fc=X(2,1)
    fa=(q*l)-fb-fc

  END SUBROUTINE find_react

  SUBROUTINE find_area_I(id,d1,d2,d3,d4,area,In)
    IMPLICIT NONE

    INTEGER, INTENT(IN) :: id
    REAL, INTENT(IN) :: d1,d2,d3,d4
    REAL, INTENT(OUT) :: area,In

    !Perform checks and calculate area and I
    IF (id == 1) THEN
       area=pi*(d1**2)
       In=(pi/4.)*(d1**4)
    ELSE IF (id == 2) THEN
       area=pi*((d2**2)-(d1**2))
       In=(pi/4.)*((d2**4)-(d1**4))
    ELSE IF (id == 3) THEN
       area=d1*d2
       In=(1./12.)*d1*(d2**3)
    ELSE IF (id == 4) THEN
       area=(d3*d4) - (d1*d2)
       In=((1./12.)*d3*(d4**3))-((1./12.)*d1*(d2**3))
    ELSE IF (id==5) THEN
       area=(2*d1*d2)+(d3*d4)
       In=(2*(((1/12.)*d1*(d2**3))+((d1*d2)*((d4/2)+d2)**2))) +  ((1/12.)*d3*(d4**3))
    ELSE
       PRINT *, "INVALID SHAPE !!"
    END IF
  END SUBROUTINE find_area_I

  SUBROUTINE find_disp_moment(pb,pc,fb,fc,In,disp,moment)
    IMPLICIT NONE

    REAL, INTENT(IN) :: pb,pc,fb,fc,In
    REAL, INTENT(OUT) :: disp(:),moment(:)

    disp=0.
    moment=0.
    !Load the displacment due to support B
    disp(2:INT(pb))=disp(2:INT(pb))+(((fb*(divl(2:INT(pb))**2))*((3*pb)-divl(2:INT(pb))))/(6.*E*In))
    disp(INT(pb)+2:lsize)=disp(INT(pb)+2:lsize)+(((fb*(pb**2))*((3*divl(INT(pb)+2:lsize))-pb))/(6.*E*In))
    !Moment due to the support B
    moment(1:INT(pb))=(moment(1:INT(pb)))+(fb*(pb-divl(1:INT(pb))))

    !Load the displacement due to support C
    disp(2:INT(pc))=disp(2:INT(pc))+(((fc*(divl(2:INT(pc))**2))*((3*pc)-divl(2:INT(pc))))/(6.*E*In))
    disp(INT(pc)+2:lsize)=disp(INT(pc)+2:lsize)+(((fc*(pc**2))*((3*divl(INT(pc)+2:lsize))-pc))/(6.*E*In))
    !Moment due to the support C
    moment(1:INT(pc))=(moment(1:INT(pc)))+(fc*(pc-divl))

    !Load the displacment due to the distributed load
    disp=disp-((q*(divl**2))*((divl**2)+(6.*(l**2))-(4.*l*divl)))/(24.*E*In)
    !Moment due to the distributed load
    moment=moment-(q*((l-divl)**2))/2.

    !Zeroing out the disp at the supports
    disp(1)=0
    disp(INT(pb)+1)=0
    disp(INT(pc)+1)=0
    moment(1)=0.
  END SUBROUTINE find_disp_moment

  SUBROUTINE find_maximum_disp_moment(disp,moment,maximum_values_disp,maximum_values_moment,maximum_disp,maximum_moment)
    IMPLICIT NONE

    REAL, INTENT(IN) :: disp(:),moment(:)
    REAL, INTENT(OUT) :: maximum_values_disp(:),maximum_values_moment(:),maximum_disp,maximum_moment

    !Finding max displacement
    maximum_values_disp(1)=ABS(MAXVAL(disp,1))
    maximum_values_disp(2)=ABS(MINVAL(disp,1))
    maximum_disp=MAXVAL(maximum_values_disp,1)

    !Finding the max moment
    maximum_values_moment(1)=ABS(MAXVAL(moment,1))
    maximum_values_moment(2)=ABS(MINVAL(moment,1))
    maximum_moment=MAXVAL(maximum_values_moment,1)

  END SUBROUTINE find_maximum_disp_moment

  !Pure function to calculate the displacement at a single point
  PURE REAL FUNCTION find_disp_point(pb,pc,fb,fc,In,point) RESULT(disp_point)
    IMPLICIT NONE

    REAL,INTENT(IN) :: pb,pc,fb,fc,In
    INTEGER,INTENT(IN) :: point

    disp_point=0.

    IF (point<=INT(pb)) THEN
       disp_point=disp_point+((fb*(divl(point)**2))*((3*pb)-divl(point)))/(6*E*In)
    ELSE IF (point>INT(pb)) THEN
       disp_point=disp_point+((fb*(pb**2))*((3*divl(point))-pb))/(6*E*In)
    END IF

    IF (point<=INT(pc)) THEN
       disp_point=disp_point+((fc*(divl(point)**2))*((3*pc)-divl(point)))/(6*E*In)
    ELSE IF (point> INT(pc)) THEN
       disp_point=disp_point+((fc*(pc**2))*((3*divl(point))-pc))/(6*E*In)
    END IF
    disp_point=disp_point-((q*(divl(point)**2))*((divl(point)**2)+(6*(l**2))-(4*l*divl(point))))/(24.*E*In)

  END FUNCTION find_disp_point


END PROGRAM bmsolve_optimized_coarrays

The error only shows up when I assign a value to the best_array locally on the coarray
Here is the error :

                                                                        
 Initial reaction at A:                     10857.156     N             
 Initial reaction at B:                     4571.4375     N             
 Initial reaction at C:                     32571.406     N             
                                                                        
 Location of maximum displacement:                              6000 mm 
                                                                        
double free or corruption (out)                                         
free(): invalid pointer                                                 
                                                                        
Program received signal SIGABRT: Process abort signal.                  
                                                                        
Backtrace for this error:                                               
                                                                        
Program received signal SIGABRT: Process abort signal.                  
                                                                        
Backtrace for this error:                                               
#0  0x7f1c13756a12 in ???                                               
#1  0x7f1c13755ba5 in ???                                               
#2  0x7f1c1356fa6f in ???                                               
#3  0x7f1c135bfc4c in ???                                               
#4  0x7f1c1356f9c5 in ???                                               
#5  0x7f1c135597f3 in ???                                               
#6  0x7f1c135b3d9d in ???                                               
#7  0x7f1c135c995b in ???                                               
#8  0x7f1c135cb79b in ???                                               
#9  0x7f1c135ce102 in ???                                               
#10  0x416d2d in ???                                                    
#11  0x4026be in ???                                                    
#12  0x7f1c1355a54f in ???                                              
#13  0x7f1c1355a608 in ???                                              
#14  0x4026f4 in ???                                                    
#15  0xffffffffffffffff in ???                                          
#0  0x7fbdedb56a12 in ???                                               
#1  0x7fbdedb55ba5 in ???                                               
#2  0x7fbded96fa6f in ???                                               
#3  0x7fbded9bfc4c in ???                                               
#4  0x7fbded96f9c5 in ???                                               
#5  0x7fbded9597f3 in ???                                               
#6  0x7fbded9b3d9d in ???                                               
#7  0x7fbded9c995b in ???                                               
#8  0x7fbded9cba7f in ???                                               
#9  0x7fbded9ce102 in ???                                               
#10  0x416d2d in ???                                                    
#11  0x4026be in ???                                                    
#12  0x7fbded95a54f in ???                                              
#13  0x7fbded95a608 in ???                                              
#14  0x4026f4 in ???                                                    
#15  0xffffffffffffffff in ???                                          
                                                                        
========================================================================
=   BAD TERMINATION OF ONE OF YOUR APPLICATION PROCESSES                
=   RANK 0 PID 59835 RUNNING AT fedora                                  
=   KILLED BY SIGNAL: 6 (Aborted)                                       
========================================================================
                                                                        
========================================================================
=   BAD TERMINATION OF ONE OF YOUR APPLICATION PROCESSES                
=   RANK 1 PID 59836 RUNNING AT fedora                                  
=   KILLED BY SIGNAL: 9 (Killed)                                        
========================================================================
Error: Command:                                                         
   `/opt/intel/oneapi/mpi/2021.6.0/bin/mpiexec -n 2 ./question_2_optimiz
failed to run.

Thank you for any help

Update: This problem is unique to Gfortran and OpenCoarrays as it doesn’t occur on the ifort compiler. I am gonna continue with ifort for now, but if anyone knows how to compile this with the CAF compiler let me know.
Update 2: So this error occurs as I’ve found due to certain pointers not being implemented yet in Gfortran and OpenCoarrays as per the Fortran wiki
https://fortranwiki.org/fortran/show/Fortran+2008+status

Single image support since 4.6. Multi-image support using OpenCoarrays (including the Fortran 2015 collective subroutines) since 5.1, except allocatable or pointer components of derived type coarrays.

1 Like
2

Here is the serial code:


PROGRAM bmsolve_optimized
  USE DISLIN
  IMPLICIT NONE

  !Compiler options: -march=skylake -mtune=skylake -O3 -mfpmath=sse -funroll-loops -fexternal -blas -fstack-arrays-lblas -fcoarray=single
  !gfortran -ldislin -I /usr/local/dislin/g95 -march=skylake -mtune=skylake -O3 -mfpmath=sse -funroll-loops -fexternal-blas -fstack-arrays -lblas -fcoarray=single /usr/local/dislin/g95/dislin.f90 question_2_optimized.f90 -o question_2_optimized_gfort
  !Possible Parrallelized version
  !gfortran -ldislin -I /usr/local/dislin/g95 -march=skylake -mtune=skylake -O3 -mfpmath=sse -funroll-loops -fexternal-blas -fstack-arrays -lblas -fcoarray=single /usr/local/dislin/g95/dislin.f90 question_2_optimized.f90 -o question_2_optimized_gfort -fopenmp

 !Data type to store details of the beam cross-section
  TYPE data
     REAL, ALLOCATABLE :: disp(:),moment(:)
     REAL :: maximum_disp,maximum_values_disp(2)
     REAL :: pb,pc
     REAL :: In,A,d1,d2,d3,d4
     REAL :: fa,fb,fc
     REAL :: maximum_normal_stress,maximum_moment,maximum_values_moment(2)
     INTEGER :: id
  END TYPE data

  !Parameters
  REAL, PARAMETER :: q=8.,disp_all=5,pi=4.0*ATAN(1.0),l=6000,E=2E5
  INTEGER,PARAMETER :: lsize=6001,dmax=500

  !Variables
  INTEGER :: i,j,k,r,t,maximum_disp_loc
  INTEGER(KIND=8)::total
  REAL, ALLOCATABLE :: div(:),divl(:)
  REAL :: t1,t2
  TYPE(data) :: s
  TYPE(data) :: best_minimum_disp,best_minimum_bends,best_arr(5),best

  !File name
  CHARACTER(LEN=*),PARAMETER :: plt_file="question_2_best_plot.gp"

  !Starting CPU clock and initializing calculation counter
  CALL CPU_TIME(t1)
  total=0

  !Allocating arrays
  ALLOCATE(div(dmax))
  ALLOCATE(divl(lsize))
  ALLOCATE(s%disp(lsize))
  ALLOCATE(s%moment(lsize))

  !Output format
  100 FORMAT (X,A,25X,A)
  110 FORMAT (X,A,G23.8,3X,A)
  120 FORMAT (X,A,G36.8,3X,A)
  130 FORMAT (X,A,G25.8,2X,A)
  140 FORMAT (X,A,G27.8,X,A)
  150 FORMAT (X,A,G26.8,2X,A)
  160 FORMAT (X,A,G30.8,X,A)
  170 FORMAT (X,A,G26.8,2X,A)
  180 FORMAT (X,A,G18.8,1X,A)
  190 FORMAT (X,A,G20.8,X,A)
  191 FORMAT (X,A,G20.8,X,A)
  192 FORMAT (X,A,G33.8,X,A)

  DO i=1,lsize
     divl(i)=REAL(i)
  END DO

  DO i=1,dmax
     div(i)=REAL(i)
  END DO


  s%pb=2000.
  s%pc=4000.

  CALL find_react(s%fa, s%fb, s%fc, s%pb, s%pc)

  PRINT *, " "
  PRINT 192, "Initial reaction at A: ",s%fa,"N"
  PRINT 192, "Initial reaction at B: ",s%fb,"N"
  PRINT 192, "Initial reaction at C: ",s%fc,"N"
  PRINT *, " "

  s%id=1
  s%d1=div(80)

  CALL find_area_I(s%id, s%d1, s%d2, s%d3, s%d4, s%A, s%In)
  CALL find_disp_moment(s%pb, s%pc, s%fb, s%fc, s%In, s%disp, s%moment)
  CALL find_maximum_disp_moment(s%disp, s%moment, s%maximum_values_disp, s%maximum_values_moment, s%maximum_disp,&
       s%maximum_moment)

  IF (MAXLOC(s%maximum_values_disp,1)==1) THEN
     maximum_disp_loc=MAXLOC(s%disp,1)
  ELSE IF (MAXLOC(s%maximum_values_disp,1)==2) THEN
     maximum_disp_loc=MINLOC(s%disp,1)
  END IF

  PRINT 192, "Location of maximum displacement: ",(maximum_disp_loc-1),"mm"
  PRINT *, " "

  !Solid circular cross-section
  s%id=1
  best_arr(1)%maximum_disp=0.
  s%disp=0.
  s%d2=0.
  s%d3=0.
  s%d4=0.

  DO i=20,dmax
     s%d1=div(i) !Radius

     CALL find_area_I(s%id, s%d1, s%d2, s%d3, s%d4, s%A, s%In)
     s%maximum_disp=ABS(find_disp_point(s%pb, s%pc, s%fb, s%fc, s%In, maximum_disp_loc))
     s%maximum_normal_stress=(s%maximum_moment*s%d1)/s%In


     IF (s%maximum_disp>best_arr(1)%maximum_disp .AND. s%maximum_disp<disp_all) THEN
        best_arr(1)=s
     END IF
     total=total+1
  END DO

  CALL find_disp_moment(best_arr(1)%pb, best_arr(1)%pc, best_arr(1)%fb, best_arr(1)%fc, best_arr(1)%In, best_arr(1)%disp,&
       best_arr(1)%moment)

  PRINT 100, "SHAPE: ","Circle"
  PRINT 110, "Moment of Inertia: ",best_arr(1)%In, "mm^4"
  PRINT 120, "Area: ",best_arr(1)%A,"mm^2"
  PRINT 130, "Maximum deflection: ",best_arr(1)%maximum_disp,"mm"
  PRINT 140, "Global maxima: ",MAXVAL(best_arr(1)%disp,1),"mm"
  PRINT 150, "Location: ",MAXLOC(best_arr(1)%disp,1)-1,"mm"
  PRINT 160, "Global minima: ", MINVAL(best_arr(1)%disp,1),"mm"
  PRINT 170, "Location: ",MINLOC(best_arr(1)%disp,1)-1,"mm"
  PRINT 180, "Maximum normal stress: ",best_arr(1)%maximum_normal_stress,"Pa"
  PRINT 190, "Global maxima moment: ",MAXVAL(best_arr(1)%moment),"N-mm"
  PRINT 191, "Global minima moment: ",MINVAL(best_arr(1)%moment),"N-mm"
  PRINT 192, "Radius: ",best_arr(1)%d1,"mm"
  PRINT *, " "
  PRINT *, "Total calculations so far: ",total
  PRINT *, " "

  s%id=2
  best_arr(2)%maximum_disp=0.
  s%disp=0.
  s%d3=0.
  s%d4=0.

  DO i=20,dmax-1
     DO j=i+1,dmax
        s%d1=div(i)
        s%d2=div(j)
        IF ((s%d2-s%d1)<10.) THEN
           CYCLE
        END IF

        CALL find_area_I(s%id, s%d1, s%d2, s%d3, s%d4, s%A, s%In)
        s%maximum_disp=ABS(find_disp_point(s%pb, s%pc, s%fb, s%fc, s%In, maximum_disp_loc))
        s%maximum_normal_stress=(s%maximum_moment*s%d2)/s%In
        IF (s%maximum_disp>best_arr(2)%maximum_disp .AND. s%maximum_disp<disp_all) THEN
           best_arr(2)=s
        END IF
        total=total+1
     END DO
  END DO

  CALL find_disp_moment(best_arr(2)%pb, best_arr(2)%pc, best_arr(2)%fb, best_arr(2)%fc, best_arr(2)%In, best_arr(2)%disp,&
       best_arr(2)%moment)

  PRINT 100, "SHAPE: ","Hollow Circle"
  PRINT 110, "Moment of Inertia: ",best_arr(2)%In, "mm^4"
  PRINT 120, "Area: ",best_arr(2)%A,"mm^2"
  PRINT 130, "Maximum deflection: ",best_arr(2)%maximum_disp,"mm"
  PRINT 140, "Global maxima: ",MAXVAL(best_arr(2)%disp,1),"mm"
  PRINT 150, "Location: ",MAXLOC(best_arr(2)%disp,1)-1,"mm"
  PRINT 160, "Global minima: ", MINVAL(best_arr(2)%disp,1),"mm"
  PRINT 170, "Location: ",MINLOC(best_arr(2)%disp,1)-1,"mm"
  PRINT 180, "Maximum normal stress: ",best_arr(2)%maximum_normal_stress,"Pa"
  PRINT 190, "Global maxima moment: ",MAXVAL(best_arr(2)%moment),"N-mm"
  PRINT 191, "Global minima moment: ",MINVAL(best_arr(2)%moment),"N-mm"
  PRINT 192, "Inner Radius: ",best_arr(2)%d1,"mm"
  PRINT 192, "Outer radius: ",best_arr(2)%d2,"mm"
  PRINT *, " "
  PRINT *, "Total calculations so far: ",total
  PRINT *, " "


  s%id=3
  best_arr(3)%maximum_disp=0.
  s%disp=0.
  s%d3=0.
  s%d4=0.


  DO i=20,dmax-1
     DO j=i+1,dmax
        s%d1=div(i)
        s%d2=div(j)

        CALL find_area_I(s%id, s%d1, s%d2, s%d3, s%d4, s%A, s%In)
        s%maximum_disp=ABS(find_disp_point(s%pb, s%pc, s%fb, s%fc, s%In, maximum_disp_loc))
        s%maximum_normal_stress=(s%maximum_moment*(s%d2/2.))/s%In

        IF (s%maximum_disp>best_arr(3)%maximum_disp .AND. s%maximum_disp<disp_all) THEN
           best_arr(3)=s
        END IF
        total=total+1
     END DO
  END DO

  CALL find_disp_moment(best_arr(3)%pb, best_arr(3)%pc, best_arr(3)%fb, best_arr(3)%fc,best_arr(3)%In, best_arr(3)%disp,&
       best_arr(3)%moment)

  PRINT 100, "SHAPE: ","Solid Rectangle"
  PRINT 110, "Moment of Inertia: ",best_arr(3)%In, "mm^4"
  PRINT 120, "Area: ",best_arr(3)%A,"mm^2"
  PRINT 130, "Maximum deflection: ",best_arr(3)%maximum_disp,"mm"
  PRINT 140, "Global maxima: ",MAXVAL(best_arr(3)%disp,1),"mm"
  PRINT 150, "Location: ",MAXLOC(best_arr(3)%disp,1)-1,"mm"
  PRINT 160, "Global minima: ", MINVAL(best_arr(3)%disp,1),"mm"
  PRINT 170, "Location: ",MINLOC(best_arr(3)%disp,1)-1,"mm"
  PRINT 180, "Maximum normal stress: ",best_arr(3)%maximum_normal_stress,"Pa"
  PRINT 190, "Global maxima moment: ",MAXVAL(best_arr(3)%moment),"N-mm"
  PRINT 191, "Global minima moment: ",MINVAL(best_arr(3)%moment),"N-mm"
  PRINT 192, "Breadth: ",best_arr(3)%d1,"mm"
  PRINT 192, "Height: ",best_arr(3)%d2,"mm"
  PRINT *, " "
  PRINT *, "Total calculations so far: ",total
  PRINT *, " "


  s%id=4
  best_arr(4)%maximum_disp=0.
  s%disp=0.

  DO i=100,dmax-1
     DO j=100,dmax-1
        DO k=i+1,dmax
           DO r=j+1,dmax

              s%d1=div(i)
              s%d2=div(j)

              s%d3=div(k)
              s%d4=div(r)

              CALL find_area_I(s%id, s%d1, s%d2, s%d3, s%d4, s%A, s%In)
              s%maximum_disp=ABS(find_disp_point(s%pb, s%pc, s%fb, s%fc, s%In, maximum_disp_loc))
              s%maximum_normal_stress=(s%maximum_moment*(s%d4/2.))/s%In

              IF (s%maximum_disp>best_arr(4)%maximum_disp .AND. s%maximum_disp<disp_all) THEN
                 best_arr(4)=s
              END IF
              total=total+1
           END DO
        END DO
     END DO
  END DO

  CALL find_disp_moment(best_arr(4)%pb, best_arr(4)%pc, best_arr(4)%fb, best_arr(4)%fc,best_arr(4)%In, best_arr(4)%disp,&
       best_arr(4)%moment)


  PRINT 100, "SHAPE: ","Hollow Rectangle"
  PRINT 110, "Moment of Inertia: ",best_arr(4)%In, "mm^4"
  PRINT 120, "Area: ",best_arr(4)%A,"mm^2"
  PRINT 130, "Maximum deflection: ",best_arr(4)%maximum_disp,"mm"
  PRINT 140, "Global maxima: ",MAXVAL(best_arr(4)%disp,1),"mm"
  PRINT 150, "Location: ",MAXLOC(best_arr(4)%disp,1)-1,"mm"
  PRINT 160, "Global minima: ", MINVAL(best_arr(4)%disp,1),"mm"
  PRINT 170, "Location: ",MINLOC(best_arr(4)%disp,1)-1,"mm"
  PRINT 180, "Maximum normal stress: ",best_arr(4)%maximum_normal_stress,"Pa"
  PRINT 190, "Global maxima moment: ",MAXVAL(best_arr(4)%moment),"N-mm"
  PRINT 191, "Global minima moment: ",MINVAL(best_arr(4)%moment),"N-mm"
  PRINT 192, "Inner Breadth: ",best_arr(4)%d1,"mm"
  PRINT 192, "Inner Height: ",best_arr(4)%d2,"mm"
  PRINT 192, "Outer Breadth: ",best_arr(4)%d3,"mm"
  PRINT 192, "Outer Height: ",best_arr(4)%d4,"mm"
  PRINT *, " "
  PRINT *, "Total calculations so far: ",total
  PRINT *, " "



  s%id=5
  best_arr(5)%maximum_disp=0.
  s%disp=0.

  DO i=10,dmax
     DO j=10,dmax
        DO k=10,dmax
           DO r=10,dmax
              s%d1=div(i)
              s%d2=div(j)

              s%d3=div(k)
              s%d4=div(r)

             IF (s%d3>s%d1 .OR. (s%d2 .NE. s%d3)) THEN
                CYCLE
             END IF

              CALL find_area_I(s%id, s%d1, s%d2, s%d3, s%d4, s%A, s%In)
              s%maximum_disp=ABS(find_disp_point(s%pb, s%pc, s%fb, s%fc, s%In, maximum_disp_loc))
              s%maximum_normal_stress=(s%maximum_moment*((s%d4+(2*s%d2))/2.))/s%In

              IF (s%maximum_disp>best_arr(5)%maximum_disp .AND. s%maximum_disp<disp_all) THEN
                 best_arr(5)=s
              END IF
              total=total+1

           END DO
        END DO
     END DO
  END DO

  CALL find_disp_moment(best_arr(5)%pb, best_arr(5)%pc, best_arr(5)%fb, best_arr(5)%fc,best_arr(5)%In, best_arr(5)%disp,&
       best_arr(5)%moment)



  PRINT 100, "SHAPE: ","I-beam"
  PRINT 110, "Moment of Inertia: ",best_arr(5)%In, "mm^4"
  PRINT 120, "Area: ",best_arr(5)%A,"mm^2"
  PRINT 130, "Maximum deflection: ",best_arr(5)%maximum_disp,"mm"
  PRINT 140, "Global maxima: ",MAXVAL(best_arr(5)%disp,1),"mm"
  PRINT 150, "Location: ",MAXLOC(best_arr(5)%disp,1)-1,"mm"
  PRINT 160, "Global minima: ", MINVAL(best_arr(5)%disp,1),"mm"
  PRINT 170, "Location: ",MINLOC(best_arr(5)%disp,1)-1,"mm"
  PRINT 180, "Maximum normal stress: ",best_arr(5)%maximum_normal_stress,"Pa"
  PRINT 190, "Global maxima moment: ",MAXVAL(best_arr(5)%moment),"N-mm"
  PRINT 191, "Global minima moment: ",MINVAL(best_arr(5)%moment),"N-mm"
  PRINT 192, "Flange Breadth: ",best_arr(5)%d1,"mm"
  PRINT 192, "Flange Height: ",best_arr(5)%d2,"mm"
  PRINT 192, "Web Breadth: ",best_arr(5)%d3,"mm"
  PRINT 192, "Web Height: ",best_arr(5)%d4,"mm"
  PRINT *, " "
  PRINT *, "Total calculations so far: ",total
  PRINT *, " "


  best%maximum_disp=0.

  DO i=1,5
     IF (best_arr(i)%maximum_disp>best%maximum_disp) THEN
        best=best_arr(i)
     ELSE IF (best_arr(i)%maximum_disp == best%maximum_disp .AND. best_arr(i)%maximum_normal_stress<best%maximum_normal_stress) THEN
        best=best_arr(i)
     END IF
  END DO

  IF (best%id==1) THEN
     PRINT *, "The best shape is the solid circle shape"
  ELSE IF (best%id==2) THEN
     PRINT *, "The best shape is the hollow circle shape"
  ELSE IF (best%id==3) THEN
     PRINT *, "The best shape is the solid rectangle shape"
  ELSE IF (best%id==4) THEN
     PRINT *, "The best shape is the hollow rectangle shape"
  ELSE IF (best%id==5) THEN
     PRINT *, "The best shape is the I-beam shape"
  END IF
  PRINT *, " "

  !Assigns the best shape to the dummy variable
  s=best
  best_minimum_disp=best
  best_minimum_bends=best


  !Finding the lowest maximum displacement by finding best place for supports
  DO i=2,lsize-1
     DO j=i+1,lsize
        s%pb=div(i)
        s%pc=divl(j)

        !Finding the reaction forces
        CALL find_react(s%fa, s%fb, s%fc, s%pb, s%pc)

        IF (s%fb .NE. s%fb .OR. s%fc .NE. s%fc) THEN
           CYCLE
        END IF

        CALL find_disp_moment(s%pb, s%pc, s%fb, s%fc, s%In, s%disp, s%moment)
        CALL find_maximum_disp_moment(s%disp, s%moment, s%maximum_values_disp, s%maximum_values_moment,&
             s%maximum_disp, s%maximum_moment)

        IF (MAXLOC(s%maximum_values_disp,1)==1) THEN
           maximum_disp_loc=MAXLOC(s%disp,1)
        ELSE IF (MAXLOC(s%maximum_values_disp,1)==2) THEN
           maximum_disp_loc=MINLOC(s%disp,1)
        END IF

        s%maximum_disp=ABS(find_disp_point(s%pb, s%pc, s%fb, s%fc, s%In, maximum_disp_loc))

        !Finding the maximum normal stress
        IF (best%id==1) THEN
           s%maximum_normal_stress=(s%maximum_moment*s%d1)/s%In
        ELSE IF (best%id==2) THEN
           s%maximum_normal_stress=(s%maximum_moment*s%d2)/s%In
        ELSE IF (best%id==3) THEN
           s%maximum_normal_stress=(s%maximum_moment*(s%d2/2.))/s%In
        ELSE IF (best%id==4) THEN
           s%maximum_normal_stress=(s%maximum_moment*(s%d4/2.))/s%In
        ELSE IF (best%id==5) THEN
           s%maximum_normal_stress=(s%maximum_moment*((s%d4+(2.*s%d2))/2.))/s%In
        END IF

        IF (s%maximum_disp<best_minimum_disp%maximum_disp) THEN
           best_minimum_disp=s
        END IF
        IF (s%maximum_normal_stress<best_minimum_bends%maximum_normal_stress) THEN
           best_minimum_bends=s
        END IF
        total=total+1
     END DO
  END DO

  PRINT *, " "
  PRINT *, "The best positions with the least displacment are: "
  PRINT 192, "Support B at: ",best_minimum_disp%pb,"mm"
  PRINT 192, "With reaction force: ",best_minimum_disp%fb,"N"
  PRINT 192, "Support C at: ",best_minimum_disp%pc,"mm"
  PRINT 192, "With reaction force: ",best_minimum_disp%fc,"N"
  PRINT 192, "Maximum displacement: ",best_minimum_disp%maximum_disp,"mm"
  PRINT 192, "Maximum normal stress: ",best_minimum_disp%maximum_normal_stress,"Pa"
  PRINT *, " "

  PRINT *, " "
  PRINT *, "The best positions with the least maximum stress are: "
  PRINT 192, "Support B at: ",best_minimum_bends%pb,"mm"
  PRINT 192, "With reaction force: ",best_minimum_bends%fb,"N"
  PRINT 192, "Support C at: ",best_minimum_bends%pc,"mm"
  PRINT 192, "With reaction force: ",best_minimum_bends%fc,"N"
  PRINT 192, "Maximum displacement: ",best_minimum_bends%maximum_disp,"mm"
  PRINT 192, "Maximum normal stress: ",best_minimum_bends%maximum_normal_stress,"Pa"
  PRINT *, " "

  PRINT *, " "
  PRINT *, "Total number of calculations are: ",total
  PRINT *, " "
  CALL CPU_TIME(t2)
  PRINT 192, "Elapsed time: ",t2-t1,"Secs"
  PRINT *, " "
  PRINT *, "Printing graphs...."

  CALL plot_with_dislin(divl, best%disp, lsize, 1,best%pb,best%pc)
  CALL plot_with_dislin(divl, best_minimum_disp%disp, lsize, 2,best_minimum_disp%pb,best_minimum_disp%pc)
  CALL plot_with_dislin(divl, best_minimum_bends%disp, lsize, 3,best_minimum_bends%pb,best_minimum_bends%pc)

  PRINT *, "Done!"

  DEALLOCATE(div); DEALLOCATE(divl); DEALLOCATE(s%disp); DEALLOCATE(s%moment)

CONTAINS

  SUBROUTINE find_react(fa,fb,fc,pb,pc)
    IMPLICIT NONE

    !Data from main program
    REAL,INTENT(IN) :: pb,pc
    REAL,INTENT(OUT) :: fa,fb,fc

    !Local variables
    REAL :: detA
    REAL :: A(2,2),A_in(2,2),X(2,1),B(2,1)


    A(1,1)=((pb**2)*((3*pb)-pb))/6.
    A(1,2)=((pb**2)*((3*pc)-pb))/6.
    A(2,1)=((pb**2)*((3*pc)-pb))/6.
    A(2,2)=((pc**2)*((3*pc)-pc))/6.

    B(1,1)=(q*(pb**2)*((6*(l**2))-(4*l*pb)+(pb**2)))/24.
    B(2,1)=(q*(pc**2)*((6*(l**2))-(4*l*pc)+(pc**2)))/24.

    detA=1/((A(1,1)*A(2,2))-(A(1,2)*A(2,1)))

    A_in(1,1)=+detA*A(2,2)
    A_in(2,1)=-detA*A(2,1)
    A_in(1,2)=-detA*A(1,2)
    A_in(2,2)=+detA*A(1,1)

    X(1,1)=(B(1,1)*A_in(1,1))+(A_in(1,2)*B(2,1))
    X(2,1)=(B(1,1)*A_in(2,1))+(A_in(2,2)*B(2,1))


    fb=X(1,1)
    fc=X(2,1)
    fa=(q*l)-fb-fc

  END SUBROUTINE find_react

  SUBROUTINE find_area_I(id,d1,d2,d3,d4,area,In)
    IMPLICIT NONE

    INTEGER, INTENT(IN) :: id
    REAL, INTENT(IN) :: d1,d2,d3,d4
    REAL, INTENT(OUT) :: area,In

    !Perform checks and calculate area and I
    IF (id == 1) THEN
       area=pi*(d1**2)
       In=(pi/4.)*(d1**4)
    ELSE IF (id == 2) THEN
       area=pi*((d2**2)-(d1**2))
       In=(pi/4.)*((d2**4)-(d1**4))
    ELSE IF (id == 3) THEN
       area=d1*d2
       In=(1./12.)*d1*(d2**3)
    ELSE IF (id == 4) THEN
       area=(d3*d4) - (d1*d2)
       In=((1./12.)*d3*(d4**3))-((1./12.)*d1*(d2**3))
    ELSE IF (id==5) THEN
       area=(2*d1*d2)+(d3*d4)
       In=(2*(((1/12.)*d1*(d2**3))+((d1*d2)*((d4/2)+d2)**2))) +  ((1/12.)*d3*(d4**3))
    ELSE
       PRINT *, "INVALID SHAPE !!"
    END IF
  END SUBROUTINE find_area_I

  SUBROUTINE find_disp_moment(pb,pc,fb,fc,In,disp,moment)
    IMPLICIT NONE

    REAL, INTENT(IN) :: pb,pc,fb,fc,In
    REAL, INTENT(OUT) :: disp(:),moment(:)

    disp=0.
    moment=0.
    !Load the displacment due to support B
    disp(2:INT(pb))=disp(2:INT(pb))+(((fb*(divl(2:INT(pb))**2))*((3*pb)-divl(2:INT(pb))))/(6.*E*In))
    disp(INT(pb)+2:lsize)=disp(INT(pb)+2:lsize)+(((fb*(pb**2))*((3*divl(INT(pb)+2:lsize))-pb))/(6.*E*In))
    !Moment due to the support B
    moment(1:INT(pb))=(moment(1:INT(pb)))+(fb*(pb-divl(1:INT(pb))))

    !Load the displacement due to support C
    disp(2:INT(pc))=disp(2:INT(pc))+(((fc*(divl(2:INT(pc))**2))*((3*pc)-divl(2:INT(pc))))/(6.*E*In))
    disp(INT(pc)+2:lsize)=disp(INT(pc)+2:lsize)+(((fc*(pc**2))*((3*divl(INT(pc)+2:lsize))-pc))/(6.*E*In))
    !Moment due to the support C
    moment(1:INT(pc))=(moment(1:INT(pc)))+(fc*(pc-divl))

    !Load the displacment due to the distributed load
    disp=disp-((q*(divl**2))*((divl**2)+(6.*(l**2))-(4.*l*divl)))/(24.*E*In)
    !Moment due to the distributed load
    moment=moment-(q*((l-divl)**2))/2.

    !Zeroing out the disp at the supports
    disp(1)=0
    disp(INT(pb)+1)=0
    disp(INT(pc)+1)=0
    moment(1)=0.
  END SUBROUTINE find_disp_moment

  SUBROUTINE find_maximum_disp_moment(disp,moment,maximum_values_disp,maximum_values_moment,maximum_disp,maximum_moment)
    IMPLICIT NONE

    REAL, INTENT(IN) :: disp(:),moment(:)
    REAL, INTENT(OUT) :: maximum_values_disp(:),maximum_values_moment(:),maximum_disp,maximum_moment

    !Finding max displacement
    maximum_values_disp(1)=ABS(MAXVAL(disp,1))
    maximum_values_disp(2)=ABS(MINVAL(disp,1))
    maximum_disp=MAXVAL(maximum_values_disp,1)

    !Finding the max moment
    maximum_values_moment(1)=ABS(MAXVAL(moment,1))
    maximum_values_moment(2)=ABS(MINVAL(moment,1))
    maximum_moment=MAXVAL(maximum_values_moment,1)

  END SUBROUTINE find_maximum_disp_moment

  !Pure function to calculate the displacement at a single point
  PURE REAL FUNCTION find_disp_point(pb,pc,fb,fc,In,point) RESULT(disp_point)
    IMPLICIT NONE

    REAL,INTENT(IN) :: pb,pc,fb,fc,In
    INTEGER,INTENT(IN) :: point

    disp_point=0.

    IF (point<=INT(pb)) THEN
       disp_point=disp_point+((fb*(divl(point)**2))*((3*pb)-divl(point)))/(6*E*In)
    ELSE IF (point>INT(pb)) THEN
       disp_point=disp_point+((fb*(pb**2))*((3*divl(point))-pb))/(6*E*In)
    END IF

    IF (point<=INT(pc)) THEN
       disp_point=disp_point+((fc*(divl(point)**2))*((3*pc)-divl(point)))/(6*E*In)
    ELSE IF (point> INT(pc)) THEN
       disp_point=disp_point+((fc*(pc**2))*((3*divl(point))-pc))/(6*E*In)
    END IF
    disp_point=disp_point-((q*(divl(point)**2))*((divl(point)**2)+(6*(l**2))-(4*l*divl(point))))/(24.*E*In)

  END FUNCTION find_disp_point

  SUBROUTINE plot_with_dislin(x,y,n,check,p_b,p_c)
    IMPLICIT NONE

    INTEGER,PARAMETER :: wp=selected_real_kind(15)
    REAL(KIND(wp)),ALLOCATABLE :: ox(:),oy(:),bx(:),cx(:),py(:)
    INTEGER,INTENT(IN) :: n, check
    REAL,INTENT(IN) :: p_b,p_c
    REAL(KIND(wp)),INTENT(IN) :: x(:),y(:)

    ALLOCATE(ox(n));ALLOCATE(oy(n));ALLOCATE(bx(3));ALLOCATE(cx(3));ALLOCATE(py(3))

    ox=0.;oy=0.
    bx=p_b;cx=p_c
    py=[MINVAL(y,1),0.,MAXVAL(y,1)]

    ! Initialisation
    CALL metafl("PDF")    ! Set the file format
    CALL setpag("DA3L")   ! Set the page type
    CALL scrmod('REVERS') ! Set the background
    CALL disini()         ! Initialize

    ! Set font
    CALL triplx

    ! Labels
    CALL name("Length (mm)","X")     ! Label for the x axis
    CALL name("Deflection (mm)","Y") ! Label for the y axis
    IF (check .EQ. 1) THEN
       CALL titlin("Deflection curve for best shape",1)
    ELSE IF (check .EQ. 2) THEN
       CALL titlin("Deflection curve for positions with least displacement",1)
    ELSE IF (check .EQ. 3) THEN
       CALL titlin("Deflection curve for positions with least stress",1)
    END IF

    ! Axis set up
    CALL labels('EXP','XYZ')
    CALL graf(x(1)-1,REAL(n)+1,x(1)-1,1000.,MINVAL(y,1),MAXVAL(y,1)+0.01,MINVAL(y,1),((ABS(MAXVAL(y,1))+0.01)+ABS(MINVAL(y,1)))/10.)

    !Create the title
    CALL title

    ! Plot the graph
    CALL incmrk(0)
    CALL lintyp(2)
    CALL curve(x,y,n)
    CALL incmrk(0)
    CALL lintyp(3)
    CALL curve(x,oy,n)
    CALL incmrk(0)
    CALL lintyp(3)
    CALL curve(bx,py,3)
    CALL incmrk(0)
    CALL lintyp(3)
    CALL curve(cx,py,3)

    !End dislin
    CALL disfin

    DEALLOCATE(ox);DEALLOCATE(oy);DEALLOCATE(bx);DEALLOCATE(cx);DEALLOCATE(py)

  END SUBROUTINE plot_with_dislin

END PROGRAM bmsolve_optimized
3

I have tried out a test program and it compiles and runs, so I am not sure why my original program has a hiccup

PROGRAM coarray_access_test
  IMPLICIT NONE

  TYPE test
     REAL :: c

  END type test

  REAL :: a[*]
  REAL :: b(2)[*]
  TYPE(test) :: t(2)[*]

  a = THIS_IMAGE()
  b(1) = 0.
  b(2) = 5.
  t(1)%c = 0.
  t(2)%c = 2.

  SYNC ALL()
  PRINT *, " "
  PRINT *, "This image ",THIS_IMAGE()," has ", a
  PRINT *, "and ", b(1)," and", b(2)
  PRINT *, "and ",t(1)%c," and",t(2)%c
  PRINT *, " "

  SYNC ALL()

END PROGRAM coarray_access_test
4

Can you compile your original code in debug mode to trace the point of failure?
I compiled your final code with ifort, and it runs successfully.
If you believe this is a bug in OpenCoarrays, please report it to the authors on GitHub. I am sure they are eager to hear about it.

1 Like
5

Yup on ifort, I had no problems. I will compile it on debug mode and post back what I get

1 Like
6

The coarray in your example program is of a type with allocatable components. Your simpler example does not, and so does not encounter the missing feature.

2 Likes
7

I’m not so sure about that. I reported what I think may be the same underlying problem 9 months ago and it’s gotten no response. My impression is that effort on opencoarrays has largely ended, and has now shifted to its successor caffeine(?) currently under construction.

2 Likes
8

You’re mostly correct. We still try and fix critical bugs or do contract work on opencoarrays. Our primary subcontractor has limited time these days as well.

Caffeine still has a ways to go and will need some support from compilers to be usable without direct calls, but that is where the primary focus is.

3 Likes