I need some help to modify this code (simple 2D model) for three layers with fault angle = 10m; nx= 1000; nz=500

3

Hi Mohammed,

Given your previous questions,

I’d really encourage you to pursue a more modular and flexible program structure, and also to search for better ways (as in geometrical representations) to describe your models. In most of your threads the convoluted if-elseif blocks and conditions are very difficult to reason about. I have the feeling like you don’t really care about this, since these are probably one-off calculations.

I tried to demonstrate to you a different approach yesterday with the point-in-polygon test. Today I’m sharing with you an implicit line representation I once found in the lectures notes on Analysis of Binary Images by Robyn Owens from the University of Western Australia.

program fault

   implicit none

   integer, parameter :: nx = 1001, nz = 501

   real, parameter :: pi = 4.*atan(1.)
   real :: depth, theta_deg, theta, x0, z0, rho
   real :: dx, dz

   character(1), allocatable :: icode(:,:)

   ! grid spacing
   dx = 1.
   dz = 1.

   allocate(icode(1:nx,1:nz))

   ! fault depth
   depth = 100.

   ! fault angle
   theta_deg = 10.                ! degrees
   theta = theta_deg * (pi/180.)  ! radian

   ! fault cuts through this point
   x0 = 500.
   z0 = -50.

!
! the fault is described by the implicit line
!
!   x*sin(phi) - z*cos(phi) + rho = 0
!
! where phi = pi - theta. The meanings of the parameters are as follows
!   
!   rho ... distance from the origin (0,0) to the closest point
!           on the line
!   phi ... angle between the line and the x-axis
!   theta ... adjacent angle to phi
!
! a full description of the math describing the line is available at
!  https://homepages.inf.ed.ac.uk/rbf/CVonline/LOCAL_COPIES/OWENS/LECT2/node3.html#SECTION00031100000000000000
!

   ! calculate line constant from the known angle, and a point
   ! lying on the fault
   rho = z0*cos(pi - theta) - x0*sin(pi - theta)

   ! fill icode mask 
   call two_layers_with_fault( &
      ndim = [nx,nz], &
      spacing = [dx,dz], &
      icode = icode, &
      rho = rho, &
      theta = theta, &
      depth = depth)

   call save_txt('model.txt',icode)
   call save_pgm('model.pgm',icode,ngray=3)

contains


   !> Fill mask array for model with two horizontal layers and a fault
   !>
   !>   The mask values are 
   !>     - 0 or 1 (top layer), depending on the side of the fault line
   !>     - 2 (bottom layer)
   !>
   subroutine two_layers_with_fault(ndim,spacing,icode,rho,theta,depth)
      integer, intent(in) :: ndim(2)  ! (nx, nz)
      real, intent(in) :: spacing(2)  ! (dx, dz)
      character(1), intent(out) :: icode(:,:)
      real, intent(in) :: rho, theta, depth
         ! geometry parameters

      real :: xx, zz
      integer :: j, i

      do j = 1, ndim(2)
         do i = 1, ndim(1)

            ! we assume implicitly that point (i=1,j=1) is the origin (0,0),
            ! increasing j corresponds to negative z-values in 
            ! physical coordinate,

            xx = (i-1) * spacing(1)
            zz = 0.0 - (j-1) * spacing(2)

            if (abs(zz) < depth) then
               ! fault layer
               icode(i,j) = merge('0','1', &
                  is_below_fault(xx,zz,rho,theta))
            else
               ! bottom layer
               icode(i,j) = '2'
            end if
         end do
      end do

   end subroutine

   !> Check if a point is below the fault line
   !> 
   !> The line must be prescribed by the implicit form
   !>  
   !>   x * sin(pi - theta) - y * cos(pi - theta) + rho = 0
   !>
   pure function is_below_fault(xx,zz,rho,theta) result(res)
      real, intent(in) :: xx, zz  ! point of interest
      real, intent(in) :: rho     ! line constant
      real, intent(in) :: theta   ! angle in radians

      real, parameter :: pi = 4.*atan(1.)

      logical :: res

      if (xx*sin(pi - theta) - zz*cos(pi - theta) + rho < 0) then
         res = .true.
      else
         res = .false.
      end if

   end function

   !> Save as dumb ASCII txt file (use .txt extension)
   subroutine save_txt(filename,icode)
      character(len=*), intent(in) :: filename
      character(1), intent(in) :: icode(:,:)

      integer :: i,j,unit

      open(newunit=unit,file=filename,status='unknown')
      do j = 1, nz
         write(unit,'(*(A1))') (icode(i,j), i = 1, nx)
      end do
      close(unit)

   end subroutine

   !> Save as plain Portable GrayMap (use .pgm extension)
   subroutine save_pgm(filename,icode,ngray)
      character(len=*), intent(in) :: filename
      character(1), intent(in) :: icode(:,:)
      integer, intent(in) :: ngray
         ! number of gray values

      integer :: i, j, unit

      open(newunit=unit,file=filename,status='unknown')

      ! a description of the format is available at
      !   https://en.wikipedia.org/wiki/Netpbm

      ! header part
      write(unit,'(A)') 'P2'
      write(unit,'(I0,1X,I0)') size(icode,1), size(icode,2)
      write(unit,'(I0)') ngray

      ! pixel data, gray values
      do j = 1, size(icode,2)
         write(unit,'(*(A1,:,1X))') (icode(i,j), i = 1, size(icode,1))
      end do

      close(unit)

   end subroutine

end program

Picking the right mathematical representation can make the job much easier for you.
Here’s the output I got in the second try, after adjusting a < symbol in the is_below_fault function.

model
model1001×501 5.63 KB

3 Likes