Intel® Fortran Compiler
Build applications that can scale for the future with optimized code designed for Intel® Xeon® and compatible processors.
29628 Discussions

Fortran and Google AI

JohnNichols
Honored Contributor I
‎08-21-2026 05:32 PM
232 Views

I asked GOOGLE Ai for a Fortran program, actually three, just to see if it was easy hard or just plain dumb. 

It is just plain dumb. 

This is a 3D plastic beam analysis program,  it works until you ask it to include the Pardiso solver and then it tanks.  I know how to do it with an older program just steal the Fortran code, but I was wondering about AI.  

Steve:  Can you make this work without relying on @mecej4 rather brilliant solution.   I cannot not.   I wish I could ask Jim Dempsey or mecej4 as well but they have disappeared from the forum. 

! =====================================================================
! 3D Space Frame Plastic Solver with Global Sparse Matrix Assembly,
! CSR Pointer Array Mapping, and an Intel MKL PARDISO Solver Engine
! =====================================================================
module global_fea_mod
    use mkl_pardiso
    implicit none
    private
    public :: element_type, setup_element, assemble_sparse_triplets, solve_sparse_pardiso

    ! System topology sizing limits
    integer, parameter, public :: max_nodes = 3
    integer, parameter, public :: dof_per_node = 6
    integer, parameter, public :: total_g_dof = max_nodes * dof_per_node
    integer, parameter, public :: max_triplets = 12 * 12 * 2 ! Max non-zeros across 2 elements

    type :: element_type
        integer :: id
        integer :: node_i, node_j           ! Node connectivity indices
        integer :: LM(12)                   ! Location Matrix mapping local DoFs to global DoFs
        real(8) :: length                    ! Element physical length
        real(8) :: cx, cy, cz                ! Transformation direction cosines
        real(8) :: P, My1, Mz1, My2, Mz2    ! Active internal element forces
        real(8) :: Py, My_yield             ! Cross-sectional plastic limits
        real(8) :: k_factor                  ! Tangent stiffness softening modifier
        real(8) :: ratio                     ! Yield surface interaction ratio
        logical :: is_plastic
    end type element_type

contains

    ! Initialize element parameters and build its unique LM mapping array
    subroutine setup_element(elem, id, ni, nj, xi, yi, zi, xj, yj, zj, py, my_y)
        type(element_type), intent(out) :: elem
        integer, intent(in) :: id, ni, nj
        real(8), intent(in) :: xi, yi, zi, xj, yj, zj, py, my_y
        real(8) :: dx, dy, dz
        integer :: d

        elem%id = id
        elem%node_i = ni
        elem%node_j = nj
        elem%Py = py
        elem%My_yield = my_y
        elem%P = 0.0d0; elem%My1 = 0.0d0; elem%Mz1 = 0.0d0; elem%My2 = 0.0d0; elem%Mz2 = 0.0d0
        elem%k_factor = 1.0d0; elem%ratio = 0.0d0; elem%is_plastic = .false.

        dx = xj - xi; dy = yj - yi; dz = zj - zi
        elem%length = sqrt(dx**2 + dy**2 + dz**2)
        elem%cx = dx / elem%length; elem%cy = dy / elem%length; elem%cz = dz / elem%length

        ! Populate Location Matrix (LM) destination vector fields
        do d = 1, 6
            elem%LM(d)     = (ni - 1) * dof_per_node + d  ! Node I Global DoFs
            elem%LM(d + 6) = (nj - 1) * dof_per_node + d  ! Node J Global DoFs
        end do
    end subroutine setup_element

    ! Formulate elemental Ke + Kg, and generate index coordinate tuples for global assembly
    subroutine assemble_sparse_triplets(elem, trip_count, rows, cols, values)
        type(element_type), intent(inout) :: elem
        integer, intent(inout) :: trip_count
        integer, intent(inout) :: rows(max_triplets), cols(max_triplets)
        real(8), intent(inout) :: values(max_triplets)
        
        real(8) :: Ke(12,12), Kg(12,12), Kt(12,12)
        real(8) :: EA, EI, L, P
        integer :: r, c, global_r, global_c

        L = elem%length
        P = elem%P
        EA = 2.0d5
        EI = 4.0d3

        ! Update localized cross-section material plasticity ratios
        elem%ratio = (elem%P / elem%Py)**2 + (max(abs(elem%My1), abs(elem%My2)) / elem%My_yield)**2
        if (elem%ratio >= 1.0d0) then
            elem%is_plastic = .true.
            elem%k_factor = 0.02d0  ! Post-yield tangent stiffness drop
        else
            elem%is_plastic = .false.
            elem%k_factor = 1.0d0   ! Fully elastic
        end if

        Ke = 0.0d0; Kg = 0.0d0

        ! 1. Local Elastic Tangent Stiffness Matrix (Ke)
        Ke(1,1)   =  (EA / L)  * elem%k_factor;  Ke(1,7)   = -(EA / L)  * elem%k_factor
        Ke(7,1)   = -(EA / L)  * elem%k_factor;  Ke(7,7)   =  (EA / L)  * elem%k_factor
        Ke(5,5)   =  (4.0d0 * EI / L) * elem%k_factor; Ke(5,11)  =  (2.0d0 * EI / L) * elem%k_factor
        Ke(11,5)  =  (2.0d0 * EI / L) * elem%k_factor; Ke(11,11) =  (4.0d0 * EI / L) * elem%k_factor

        ! 2. Local Geometric Stiffness matrix formulation (Kg) for axial P-Delta tracking
        Kg(5,5)   =  (2.0d0 * P * L / 15.0d0);  Kg(5,11)  = -(P * L / 30.0d0)
        Kg(11,5)  = -(P * L / 30.0d0);          Kg(11,11) =  (2.0d0 * P * L / 15.0d0)

        Kt = Ke + Kg

        ! 3. Map sparse coordinate rows/cols using Location Matrix (LM) pointers
        do r = 1, 12
            global_r = elem%LM(r)
            do c = 1, 12
                global_c = elem%LM(c)
                
                if (global_r <= total_g_dof .and. global_c <= total_g_dof) then
                    trip_count = trip_count + 1
                    rows(trip_count) = global_r
                    cols(trip_count) = global_c
                    values(trip_count) = Kt(r, c)
                end if
            end do
        end do
    end subroutine assemble_sparse_triplets

    ! Modernized Intel MKL PARDISO Engine wrapper
    subroutine solve_sparse_pardiso(n, nnz, ia, ja, a, rhs, displacement)
        integer, intent(in) :: n, nnz
        integer, intent(in) :: ia(n+1), ja(nnz)
        real(8), intent(in) :: a(nnz)
        real(8), intent(inout) :: rhs(n)
        real(8), intent(out) :: displacement(n)

        ! PARDISO control variables
        type(MKL_PARDISO_HANDLE) :: pt(64)
        integer :: maxfct, mnum, mtype, phase, nrhs, error, msglvl
        integer :: iparm(64)
        integer :: idum(1)
        real(8) :: ddum(1)
        integer :: i

        ! Initialize PARDISO internal state pointers
        do i = 1, 64
            pt(i)%DUMMY = 0
            iparm(i) = 0
        end do

        ! Set PARDISO parameters for real structurally symmetric indefinite arrays
        iparm(1) = 1       ! Explicit parameters configuration (no defaults)
        iparm(2) = 2       ! Parallel fill-in reordering from METIS
        iparm(3) = 0       ! Automatically query environment OpenMP threads
        iparm(4) = 0       ! Direct solver execution (no iterative refinement)
        iparm(6) = 0       ! Output solution directly into displacement array
        iparm(10) = 13     ! Pivot perturbation tolerance threshold
        iparm(11) = 1      ! Enable scaling/equilibration vectors
        iparm(18) = -1     ! Output number of non-zero factorization profiles
        iparm(21) = 1      ! Bunch-Kaufman pivoting for indefinite structures

        mtype = -2         ! Real structurally symmetric matrix configuration
        nrhs = 1           ! Number of right-hand side force vectors
        maxfct = 1         ! Maximum number of numerical factorizations
        mnum = 1           ! Matrix number in solver memory stack
        msglvl = 0         ! Silent diagnostic messaging flag (set to 1 to debug)

        ! Step 1: Symbolic Analysis / Reordering Phase
        phase = 11
        call pardiso(pt, maxfct, mnum, mtype, phase, n, a, ia, ja, idum, nrhs, iparm, msglvl, ddum, ddum, error)
        if (error /= 0) then
            print *, "[PARDISO ERROR] Symbolic analysis failure code: ", error
            return
        end if

        ! Step 2: Numerical Factorization Phase
        phase = 22
        call pardiso(pt, maxfct, mnum, mtype, phase, n, a, ia, ja, idum, nrhs, iparm, msglvl, ddum, ddum, error)
        if (error /= 0) then
            print *, "[PARDISO ERROR] Matrix factorization failure code: ", error
            return
        end if

        ! Step 3: Back Substitution and Solution Vector Evaluation Phase
        phase = 33
        call pardiso(pt, maxfct, mnum, mtype, phase, n, a, ia, ja, idum, nrhs, iparm, msglvl, rhs, displacement, error)
        if (error /= 0) then
            print *, "[PARDISO ERROR] Back-substitution execution failure code: ", error
            return
        end if

        ! Step 4: Solver Memory Release Cleanup Phase
        phase = -1
        call pardiso(pt, maxfct, mnum, mtype, phase, n, ddum, ia, ja, idum, nrhs, iparm, msglvl, ddum, ddum, error)

    end subroutine solve_sparse_pardiso

end module global_fea_mod

! Main structural solver execution routine
program main
    use global_fea_mod
    implicit none

    type(element_type) :: framework(2)
    integer :: csv_unit, step, i, j, r, c
    integer, parameter :: total_steps = 50

    ! Sparse triplet tracking structures
    integer :: triplet_count
    integer :: trip_rows(max_triplets), trip_cols(max_triplets)
    real(8) :: trip_vals(max_triplets)

    ! Final CSR storage matrices passed into PARDISO
    integer :: csr_ia(total_g_dof + 1)
    integer, allocatable :: csr_ja(:)
    real(8), allocatable :: csr_a(:)
    real(8) :: dense_temp(total_g_dof, total_g_dof)

    real(8) :: F_global(total_g_dof)
    real(8) :: U_global(total_g_dof)
    
    ! Displacement control tracking parameters
    real(8) :: lambda, reference_P, reference_M, target_displacement_increment
    integer :: active_nnz

    ! Node 1: (0,0,0) [Fixed base], Node 2: (4,0,0) [Joint], Node 3: (8,3,0) [Free end]
    call setup_element(framework(1), 1, 1, 2, 0.0d0, 0.0d0, 0.0d0, 4.0d0, 0.0d0, 0.0d0, 600.0d0, 90.0d0)
    call setup_element(framework(2), 2, 2, 3, 4.0d0, 0.0d0, 0.0d0, 8.0d0, 3.0d0, 0.0d0, 450.0d0, 70.0d0)

    reference_P = 20.0d0
    reference_M = 4.5d0
    target_displacement_increment = 0.005d0     ! Prescribe fixed 5mm structural movement steps
    lambda = 0.0d0                              

    open(newunit=csv_unit, file='global_assembly_history.csv', status='replace', action='write')
    write(csv_unit, '(A)') "Step,Lambda,ElemID,Length,P_Internal,My_Internal,Node2_U_Disp,Node2_V_Rot,YieldRatio,State"

    print *, "====================================================================================="
    print *, "            3D ADAPTIVE FEA FRAME ENGINE - INTEL MKL PARDISO SPARSE SOLVER           "
    print , "====================================================================================="
    write(, '(A5, A8, A4, A12, A12, A12, A8)') "Step", "Lambda", "ID", "P-Internal", "Node2-Disp", "Yield-R", "Status"
    print *, "-------------------------------------------------------------------------------------"
    ! Execution step loop
    do step = 1, total_steps
        triplet_count = 0
        F_global = 0.0d0
        dense_temp = 0.0d0
        ! 1. Assemble elements into sparse triplet coordinate vectors
        do i = 1, 2
            call assemble_sparse_triplets(framework(i), triplet_count, trip_rows, trip_cols, trip_vals)
        
            ! 2. Compress triplets onto a temporary matrix to eliminate duplicate index slots
            do i = 1, triplet_count
                dense_temp(trip_rows(i), trip_cols(i)) = dense_temp(trip_rows(i), trip_cols(i)) + trip_vals(i)
            end do
        ! 3. ADAPTIVE DISPLACEMENT CONTROL BLOCK: Evaluate active diagonal tangent stiffness
        if (dense_temp(7,7) > 1.0d-4) then
            lambda = lambda + (dense_temp(7,7) * target_displacement_increment) / reference_P
        else
            lambda = lambda + 0.01d0
        end if
        ! Map loads based on adaptive load factors
        F_global(7) = reference_P * lambda
        F_global(11) = reference_M * lambda
        ! 4. Convert structural dense components into standard Compressed Sparse Row (CSR) format
        active_nnz = 0
        do r = 1, total_g_dof
            do c = 1, total_g_dof
                if (abs(dense_temp(r,c)) > 1.0d-10) then
                    active_nnz = active_nnz + 1
                end if
            end do
        end do
    end do
    
    if (allocated(csr_ja)) deallocate(csr_ja)
    if (allocated(csr_a)) deallocate(csr_a)
    allocate(csr_ja(active_nnz), csr_a(active_nnz))
    active_nnz = 1
    do r = 1, total_g_dof
        csr_ia(r) = active_nnz
        do c = 1, total_g_dof
            if (abs(dense_temp(r,c)) > 1.0d-10) then
                csr_a(active_nnz) = dense_temp(r,c)
                csr_ja(active_nnz) = c
                active_nnz = active_nnz + 1
            end if
        end do
    end do
    csr_ia(total_g_dof + 1) = active_nnz
    ! 5. Call the multi-threaded Intel PARDISO Direct Sparse Solver module
    call solve_sparse_pardiso(total_g_dof, active_nnz - 1, csr_ia, csr_ja, csr_a, F_global, U_global)
    
    ! Post-convergence state processing
    do i = 1, 2
        framework(i)%P   = reference_P * lambda
        framework(i)%My1 = reference_M * lambda
        write(csv_unit, '(I0,A,F6.2,A,I0,A,F6.2,A,F7.1,A,F7.1,A,E11.4,A,E11.4,A,F6.2,A,A)') &step, ',', lambda, ',', framework(i)%id, ',', framework(i)%length, ',', &framework(i)%P, ',', framework(i)%My1, ',', U_global(7), ',', U_global(11), ',', &framework(i)%ratio, ',', merge("PLASTIC", "ELASTIC", framework(i)%is_plastic)
        write(*, '(I4, F8.2, I4, F12.1, E12.3, F12.2, 2X, A7)') step, lambda, framework(i)%id, framework(i)%P, &U_global(7), framework(i)%ratio, &merge("PLASTIC", "ELASTIC", framework(i)%is_plastic)
    end do
    end do
    close(csv_unit)
    print *, "-------------------------------------------------------------------------------------"
    print *, "PARDISO execution complete! Sparse matrix tracking saved to 'global_assembly_history.csv'."
    end program main

 

!****************************************************************************
!
!  PROGRAM: Sulphur
!
!  PURPOSE:  Entry point for the console application.
!
!   Technical Implementation Details
!   Passivation ModelingMechanism: 
!   As solid {CaSO}_{4}) precipitates, it adheres to the solid {CaCO}_{3} particle core.
!   Equation: passivation_factor = exp(-beta * state(3))
!   Behavior: The reaction slows down exponentially relative to the product mass. 
!   It mimics physical shell blinding.4th-Order Runge-Kutta 
!   SolverPrecision: RK4 cuts the localized truncation error down to (O(dt^5)).
!   Process: It samples the derivative curves four distinct times across each time step segment.
!   Stability: It avoids numerical oscillations or artificial negative masses when a reaction terminates abruptly.
!****************************************************************************

 program caco3_h2so4_rk4_passivation
    implicit none

    ! Molar masses (g/mol)
    real, parameter :: MM_CaCO3 = 100.09
    real, parameter :: MM_H2SO4 = 98.08
    real, parameter :: MM_CaSO4 = 136.14
    real, parameter :: MM_CO2   = 44.01

    ! Physical and Thermodynamic Constants
    real, parameter :: R = 8.314          
    real, parameter :: A_freq = 1.2e4      
    real, parameter :: Ea = 35000.0        
    real, parameter :: volume = 1.0       

    ! Passivation Parameter
    ! Higher values simulate faster surface choking by insoluble CaSO4
    real, parameter :: beta = 0.08         

    ! Simulation Settings
    real :: temperature                    
    real :: t, dt, t_max
    integer :: step, total_steps
    integer, parameter :: file_unit = 15   

    ! State Vector: [1]=CaCO3, [2]=H2SO4, [3]=CaSO4, [4]=CO2
    real :: y(4)
    real :: k1(4), k2(4), k3(4), k4(4), y_tmp(4), dydt(4)
    
    real :: k_arrhenius                    

    ! 1. Initialize State and Conditions
    temperature = 298.15   ! 25 °C
    y(1) = 5000.0            ! Initial CaCO3 (g)
    y(2) = 49.0            ! Initial H2SO4 (g)
    y(3) = 0.0             ! Initial CaSO4 (g)
    y(4) = 0.0             ! Initial CO2 (g)

    ! Time Configuration
    t = 0.0
    dt = 0.1               ! RK4 allows a larger, stable time step than Euler
    t_max = 60.0           
    total_steps = int(t_max / dt)

    ! Calculate global Arrhenius constant
    k_arrhenius = A_freq * exp(-Ea / (R * temperature))

    ! 2. Setup File Export
    open(unit=file_unit, file='passivation_kinetics.csv', status='replace', action='write')
    write(file_unit, '(A)') "Time_s,CaCO3_g,H2SO4_g,CaSO4_g,CO2_g"

    print '(A, F6.2, A)', "--> Simulating Passivation with RK4 at: ", temperature - 273.15, " C"
    print '(A10, A12, A12, A12, A12)', "Time(s)", "CaCO3(g)", "H2SO4(g)", "CaSO4(g)", "CO2(g)"
    print '(A)', "------------------------------------------------------------"

    ! 3. RK4 Integration Loop
    do step = 0, total_steps
        
        ! Export current state to CSV
        write(file_unit, '(F8.2, A1, F10.3, A1, F10.3, A1, F10.3, A1, F10.3)') &
            t, ',', y(1), ',', y(2), ',', y(3), ',', y(4)

        ! Print to console every 5 seconds
        if (mod(step, 50) == 0) then
            print '(F10.2, F12.2, F12.2, F12.2, F12.2)', t, y(1), y(2), y(3), y(4)
        end if

        ! Early exit conditions
        if (y(1) <= 0.001 .or. y(2) <= 0.001) then
            print '(A, F5.2, A)', "--> Reaction stopped at ", t, " s (reactant depletion or passivation halt)."
            exit
        end if

        ! RK4 Step Evaluation
        ! Evaluation 1
        call get_derivatives(y, dydt, k_arrhenius, beta, volume)
        k1 = dydt * dt

        ! Evaluation 2
        y_tmp = max(0.0, y + 0.5 * k1)
        call get_derivatives(y_tmp, dydt, k_arrhenius, beta, volume)
        k2 = dydt * dt

        ! Evaluation 3
        y_tmp = max(0.0, y + 0.5 * k2)
        call get_derivatives(y_tmp, dydt, k_arrhenius, beta, volume)
        k3 = dydt * dt

        ! Evaluation 4
        y_tmp = max(0.0, y + k3)
        call get_derivatives(y_tmp, dydt, k_arrhenius, beta, volume)
        k4 = dydt * dt

        ! Final Weighted Update
        y = max(0.0, y + (k1 + 2.0*k2 + 2.0*k3 + k4) / 6.0)

        t = t + dt
    end do

    close(file_unit)
    print '(A)', "------------------------------------------------------------"
    print '(A)', "--> Simulation complete. Saved to 'passivation_kinetics.csv'"

contains

    ! Subroutine to calculate derivatives for the state vector
    subroutine get_derivatives(state, rates, k_arr, b_coeff, vol)
        real, intent(in)  :: state(4)
        real, intent(out) :: rates(4)
        real, intent(in)  :: k_arr, b_coeff, vol
        real :: c_h2so4, r_molar, passivation_factor

        ! Safeguard against negative mass states during intermediate RK4 steps
        if (state(1) <= 0.0 .or. state(2) <= 0.0) then
            rates = 0.0
            return
        end if

        ! Calculate acid concentration (mol/L)
        c_h2so4 = (state(2) / MM_H2SO4) / vol

        ! Passivation Model: Exponential decay of active surface area due to solid CaSO4 build-up
        passivation_factor = exp(-b_coeff * state(3))

        ! Total Molar Rate combining Shrinking Core and Passivation Blockage
        r_molar = k_arr * (state(1)**(2.0/3.0)) * c_h2so4 * passivation_factor

        ! Convert to Mass Rates (g/s)
        rates(1) = -r_molar * MM_CaCO3  ! d(CaCO3)/dt
        rates(2) = -r_molar * MM_H2SO4  ! d(H2SO4)/dt
        rates(3) =  r_molar * MM_CaSO4  ! d(CaSO4)/dt
        rates(4) =  r_molar * MM_CO2    ! d(CO2)/dt
    end subroutine get_derivatives

end program caco3_h2so4_rk4_passivation

This one models limestone with impurities, as long as you are not over about 150 lines it is ok, once past that it is a world of hurt. 

 

0 Kudos

Link Copied

4 Replies
JohnNichols
Honored Contributor I
‎08-21-2026 05:35 PM
231 Views

AI makes stupid mistakes such as 

integer k
REAL K(30)

That one occurred in 2 of the three programs.  

When you point out the error politely, it says, well it came in the supplied code or it should have known or some such excuse.  

GOOGLE AI also suggests GFORTRAN.  

0 Kudos
Copy link
Steve_Lionel
Honored Contributor III
‎08-22-2026 08:27 AM
181 Views

This sort of thing is not my area of expertise. If the program failed to compile, that would be another thing. I am not a fan of having AI write code for you - it tends to look good but be overly complex and have numerous flaws. Often it will hallucinate functions and libraries that don't exist.

0 Kudos
Copy link
JohnNichols
Honored Contributor I
‎08-22-2026 10:21 AM
159 Views
In response to Steve_Lionel

Steve:

I agree with you, it proved to be a circular task, and yes numerous flaws.  

But I was just interested to understand what it produced.  It is fixable, but the cost would be as high as just developing existing code or writing it from scratch.  Anyone who wastes corporate funds on this is not going to get a good return.  

That is why I miss the two amigos so much.  

In response to Steve_Lionel
0 Kudos
Copy link
cean
New Contributor II
‎08-24-2026 05:25 AM
37 Views 
        do i = 1, 2
            call assemble_sparse_triplets(framework(i), triplet_count, trip_rows, trip_cols, trip_vals)
        
            ! 2. Compress triplets onto a temporary matrix to eliminate duplicate index slots
            do i = 1, triplet_count

check line 236-242.

0 Kudos
Copy link
Reply

Community support is provided Monday to Friday. Other contact methods are available here.

Intel does not verify all solutions, including but not limited to any file transfers that may appear in this community. Accordingly, Intel disclaims all express and implied warranties, including without limitation, the implied warranties of merchantability, fitness for a particular purpose, and non-infringement, as well as any warranty arising from course of performance, course of dealing, or usage in trade.