Fortran Source Code & Offline Testing

Grid Jacobian Ratio

Core Numerical Engine in Fortran 90 • 27 total downloads

jacobian.f90

! =========================================================================
! Source File: jacobian.f90
! =========================================================================

!==============================================================================
! ThermoFluidCalc — Calculator #29 : Jacobian of Transformation
!==============================================================================
! J = x_xi * y_eta - x_eta * y_xi
! Bilinear mapping (xi,eta in [0,1]):
!   x = (1-xi)(1-eta)x1 + xi(1-eta)x2 + xi*eta*x3 + (1-xi)*eta*x4
! Reference : Gupta, §5.6
! Build: gfortran -O2 -o jacobian jacobian.f90
!==============================================================================
program jacobian
  implicit none
  integer, parameter :: dp = selected_real_kind(15,307)
  integer, parameter :: MX=200
  integer :: mode,N,M,i,j,idx
  real(dp) :: vx(4),vy(4),xi,eta
  real(dp) :: x_xi,x_eta,y_xi,y_eta,Jd
  real(dp) :: Jc,Jcorner(4),Jmin,Jmax,Jrat
  real(dp) :: xi_x,xi_y,eta_x,eta_y
  real(dp) :: gx(0:MX,0:MX),gy(0:MX,0:MX)
  real(dp) :: cx,cy,Jarr(MX*MX),s1,s2,avg,std,mn,mxx
  integer :: neg
  character(len=20) :: rating
  real(dp) :: cxi(4),ceta(4)

  data cxi  /0.0_dp,1.0_dp,1.0_dp,0.0_dp/
  data ceta /0.0_dp,0.0_dp,1.0_dp,1.0_dp/

  read(*,*) mode

  select case(mode)
  !== MODE 1 ==
  case(1)
    read(*,*) vx(1),vy(1),vx(2),vy(2),vx(3),vy(3),vx(4),vy(4)
    ! Center
    call eval_jac(vx,vy,0.5_dp,0.5_dp,x_xi,x_eta,y_xi,y_eta,Jc)
    ! Corners
    do i=1,4
      call eval_jac(vx,vy,cxi(i),ceta(i),x_xi,x_eta,y_xi,y_eta,Jcorner(i))
    end do
    Jmin=min(Jc,Jcorner(1),Jcorner(2),Jcorner(3),Jcorner(4))
    Jmax=max(Jc,Jcorner(1),Jcorner(2),Jcorner(3),Jcorner(4))
    if(abs(Jmax)>1e-30_dp) then; Jrat=Jmin/Jmax; else; Jrat=0; end if
    call get_rat(Jrat,Jmin,rating)
    ! Full matrix at center
    call eval_jac(vx,vy,0.5_dp,0.5_dp,x_xi,x_eta,y_xi,y_eta,Jd)
    if(abs(Jd)>1e-30_dp) then
      xi_x= y_eta/Jd; xi_y=-x_eta/Jd; eta_x=-y_xi/Jd; eta_y= x_xi/Jd
    else
      xi_x=0;xi_y=0;eta_x=0;eta_y=0
    end if

    write(*,'(A,I1)')     'MODE=',mode
    write(*,'(A)')        'MODE_NAME=Single Quad'
    write(*,'(A,ES15.8)') 'J_CENTER=',Jc
    write(*,'(A,ES15.8)') 'J_MIN=',Jmin
    write(*,'(A,ES15.8)') 'J_MAX=',Jmax
    write(*,'(A,F12.6)')  'J_RATIO=',Jrat
    write(*,'(A,A)')      'RATING=',trim(rating)
    write(*,'(A,ES15.8)') 'X_XI=',x_xi
    write(*,'(A,ES15.8)') 'X_ETA=',x_eta
    write(*,'(A,ES15.8)') 'Y_XI=',y_xi
    write(*,'(A,ES15.8)') 'Y_ETA=',y_eta
    write(*,'(A,ES15.8)') 'XI_X=',xi_x
    write(*,'(A,ES15.8)') 'XI_Y=',xi_y
    write(*,'(A,ES15.8)') 'ETA_X=',eta_x
    write(*,'(A,ES15.8)') 'ETA_Y=',eta_y

    write(*,'(A)') 'CORNERS_START'
    write(*,'(A,ES15.8)') 'SW,',Jcorner(1)
    write(*,'(A,ES15.8)') 'SE,',Jcorner(2)
    write(*,'(A,ES15.8)') 'NE,',Jcorner(3)
    write(*,'(A,ES15.8)') 'NW,',Jcorner(4)
    write(*,'(A,ES15.8)') 'CTR,',Jc
    write(*,'(A)') 'CORNERS_END'

    write(*,'(A)') 'VERTS_START'
    do i=1,4; write(*,'(F12.6,A,F12.6)') vx(i),',',vy(i); end do
    write(*,'(A)') 'VERTS_END'

  !== MODE 2 ==
  case(2)
    backspace(5); read(*,*) mode,N,M
    if(N<1.or.N>MX.or.M<1.or.M>MX) then
      write(*,'(A)') 'ERROR=N,M must be 1-200.'; stop
    end if
    do j=0,M; do i=0,N; read(*,*) gx(i,j),gy(i,j); end do; end do

    idx=0; neg=0; s1=0; s2=0; mn=1e30_dp; mxx=-1e30_dp
    write(*,'(A,I1)')  'MODE=',mode
    write(*,'(A)')     'MODE_NAME=Grid Field'
    write(*,'(A,I5)')  'NI=',N
    write(*,'(A,I5)')  'NJ=',M
    write(*,'(A)') 'DATA_START'
    do j=0,M-1; do i=0,N-1
      vx(1)=gx(i,j);   vy(1)=gy(i,j)
      vx(2)=gx(i+1,j); vy(2)=gy(i+1,j)
      vx(3)=gx(i+1,j+1);vy(3)=gy(i+1,j+1)
      vx(4)=gx(i,j+1); vy(4)=gy(i,j+1)
      call eval_jac(vx,vy,0.5_dp,0.5_dp,x_xi,x_eta,y_xi,y_eta,Jd)
      cx=0.25_dp*(vx(1)+vx(2)+vx(3)+vx(4))
      cy=0.25_dp*(vy(1)+vy(2)+vy(3)+vy(4))
      idx=idx+1; Jarr(idx)=Jd
      if(Jd<0) neg=neg+1
      s1=s1+Jd; s2=s2+Jd**2
      if(Jd<mn) mn=Jd; if(Jd>mxx) mxx=Jd
      write(*,'(I4,A,I4,A,F12.6,A,F12.6,A,ES15.8)') i,',',j,',',cx,',',cy,',',Jd
    end do; end do
    write(*,'(A)') 'DATA_END'
    avg=s1/real(idx,dp); std=sqrt(max(0.0_dp,s2/real(idx,dp)-avg**2))
    write(*,'(A,I6)')    'NCELLS=',idx
    write(*,'(A,ES15.8)')'AVG_J=',avg
    write(*,'(A,ES15.8)')'MIN_J=',mn
    write(*,'(A,ES15.8)')'MAX_J=',mxx
    write(*,'(A,ES15.8)')'STD_J=',std
    write(*,'(A,I6)')    'NEG_COUNT=',neg

  !== MODE 3 ==
  case(3)
    read(*,*) vx(1),vy(1),vx(2),vy(2),vx(3),vy(3),vx(4),vy(4),xi,eta
    call eval_jac(vx,vy,xi,eta,x_xi,x_eta,y_xi,y_eta,Jd)
    cx=(1-xi)*(1-eta)*vx(1)+xi*(1-eta)*vx(2)+xi*eta*vx(3)+(1-xi)*eta*vx(4)
    cy=(1-xi)*(1-eta)*vy(1)+xi*(1-eta)*vy(2)+xi*eta*vy(3)+(1-xi)*eta*vy(4)
    if(abs(Jd)>1e-30_dp) then
      xi_x=y_eta/Jd; xi_y=-x_eta/Jd; eta_x=-y_xi/Jd; eta_y=x_xi/Jd
    else
      xi_x=0;xi_y=0;eta_x=0;eta_y=0
    end if
    write(*,'(A,I1)')     'MODE=',mode
    write(*,'(A)')        'MODE_NAME=Probe Point'
    write(*,'(A,F10.6)')  'XI=',xi
    write(*,'(A,F10.6)')  'ETA=',eta
    write(*,'(A,F12.6)')  'PHYS_X=',cx
    write(*,'(A,F12.6)')  'PHYS_Y=',cy
    write(*,'(A,ES15.8)') 'J=',Jd
    write(*,'(A,ES15.8)') 'X_XI=',x_xi
    write(*,'(A,ES15.8)') 'X_ETA=',x_eta
    write(*,'(A,ES15.8)') 'Y_XI=',y_xi
    write(*,'(A,ES15.8)') 'Y_ETA=',y_eta
    write(*,'(A,ES15.8)') 'XI_X=',xi_x
    write(*,'(A,ES15.8)') 'XI_Y=',xi_y
    write(*,'(A,ES15.8)') 'ETA_X=',eta_x
    write(*,'(A,ES15.8)') 'ETA_Y=',eta_y
    if(Jd>0) then; write(*,'(A)') 'VALID=YES'; else; write(*,'(A)') 'VALID=NO'; end if
    write(*,'(A)') 'VERTS_START'
    do i=1,4; write(*,'(F12.6,A,F12.6)') vx(i),',',vy(i); end do
    write(*,'(A)') 'VERTS_END'

  case default
    write(*,'(A)') 'ERROR=Invalid mode (1-3).'; stop
  end select

contains
  subroutine eval_jac(vx,vy,xi,eta,xxi,xet,yxi,yet,Jd)
    real(dp),intent(in)::vx(4),vy(4),xi,eta
    real(dp),intent(out)::xxi,xet,yxi,yet,Jd
    xxi=-(1-eta)*vx(1)+(1-eta)*vx(2)+eta*vx(3)-eta*vx(4)
    xet=-(1-xi)*vx(1)-xi*vx(2)+xi*vx(3)+(1-xi)*vx(4)
    yxi=-(1-eta)*vy(1)+(1-eta)*vy(2)+eta*vy(3)-eta*vy(4)
    yet=-(1-xi)*vy(1)-xi*vy(2)+xi*vy(3)+(1-xi)*vy(4)
    Jd=xxi*yet-xet*yxi
  end subroutine

  subroutine get_rat(ratio,jmin,r)
    real(dp),intent(in)::ratio,jmin
    character(len=20),intent(out)::r
    if(jmin<0) then; r='Folded/Invalid'
    else if(ratio>0.5_dp) then; r='Excellent'
    else if(ratio>0.3_dp) then; r='Good'
    else if(ratio>0.1_dp) then; r='Acceptable'
    else; r='Poor'
    end if
  end subroutine
end program jacobian

Solver Description

Determine the determinant Jacobian ratio of finite element and finite volume grid mesh elements.

Key Numerical Methods & Architecture

Input Redirection: Reads parameters sequentially from standard input (`stdin`) using Fortran sequential read (`read(*,*)`), ensuring modular integration.
Modular Design: Formulated using pure mathematical routines, separation of equations from output formatting, and precise numerical solvers (e.g. bisection, Newton-Raphson).
Standard Compliant: Written in clean, standards-compliant Fortran 90 to ensure cross-compiler compatibility.

🛠️ Local Compilation

To test this code on your machine, compile the source code file(s) using a standard Fortran compiler (e.g., `gfortran`).

Compilation Command:

gfortran -O3 jacobian.f90 -o jacobian

Execution Command:

Execute the program by feeding the sample input file into the program using stdin redirection:

jacobian < input.txt

📥 Downloads & Local Files

Preview of the required input file (input.txt):

! Mode (1=Single Quadrilateral, 2=Grid, 3=Single Point evaluation)\nCorner coordinates (x1 y1 x2 y2 x3 y3 x4 y4)
1
! Parameter 2
0.0 0.0 1.0 0.0 1.0 1.0 0.0 1.0

📦 Download Offline Bundle (ZIP) 8 dl 📄 Download input.txt 9 dl 📥 Download jacobian.f90 10 dl

💻 Fortran Source Code Library

Grid Jacobian Ratio

Solver Description

Key Numerical Methods & Architecture

🛠️ Local Compilation

📥 Downloads & Local Files