From ae30af19efcde3d12bf1692cf8a1e80b307ebab2 Mon Sep 17 00:00:00 2001
From: s2006749 <s2006749@ed.ac.uk>
Date: Fri, 20 Oct 2023 14:17:20 +0100
Subject: [PATCH] fixed memory allocation bug

---
 src/BC-Ellipsoid.f90 |  7 +++++--
 src/forces.f90       | 10 ++++++++++
 src/ibm.f90          |  4 ++--
 src/xcompact3d.f90   |  2 +-
 4 files changed, 18 insertions(+), 5 deletions(-)

diff --git a/src/BC-Ellipsoid.f90 b/src/BC-Ellipsoid.f90
index 740d2b77..b8626eeb 100644
--- a/src/BC-Ellipsoid.f90
+++ b/src/BC-Ellipsoid.f90
@@ -43,7 +43,7 @@ subroutine geomcomplex_ellip(epsi,nxi,nxf,ny,nyi,nyf,nzi,nzf,dx,yp,remp)
         zeromach = zeromach/two
     end do
     zeromach = ten*zeromach
-
+    is_inside=.false.
     !  orientation=[oriw, orii, orij, orik]
     call NormalizeQuaternion(orientation)
     !  shape=[shx, shy, shz]
@@ -83,11 +83,12 @@ subroutine geomcomplex_ellip(epsi,nxi,nxf,ny,nyi,nyf,nzi,nzf,dx,yp,remp)
             call is_inside_ellipsoid(point, position, orientation, shape, ra, zeromach, is_inside)
             !  r=sqrt_prec((xm-cexx)**two+(ym-ceyy)**two+(zm-cezz)**two)
             !  r=sqrt_prec((xm-cexx)**two+(ym-ceyy)**two)
-
             if (.not.is_inside) then
                 !  write(*,*) i, j, k
                 cycle
             endif
+            ! write(*,*) is_inside
+
             !  write(*,*) i, j, k, zm
             epsi(i,j,k)=remp
             !  write(*,*) remp
@@ -243,6 +244,8 @@ subroutine init_ellip (ux1,uy1,uz1,phi1)
     real(mytype) :: y,um,eqr
     integer :: k,j,i,ii,is,code
 
+    write(*,*) 'INSIDE INIT ELLIP'
+
     eqr=(shx*shy*shz)**(1.0/3.0)
     shape=[shx/eqr,shy/eqr,shz/eqr]
 
diff --git a/src/forces.f90 b/src/forces.f90
index 516bb221..c95fa957 100644
--- a/src/forces.f90
+++ b/src/forces.f90
@@ -40,6 +40,8 @@ contains
 
     integer :: iv,stp1,stp2,h
 
+    write(*,*) 'Inside INIT_FORCES'
+
     call alloc_x(ux01)
     call alloc_x(uy01)
     call alloc_x(ux11)
@@ -55,12 +57,16 @@ contains
     uz01 = zero
     uz11 = zero    
 
+    write(*,*) 'Alloc_x called'
+
     allocate(icvlf(nvol), icvrt(nvol), jcvlw(nvol), jcvup(nvol), zcvlf(nvol), zcvrt(nvol))
     allocate(icvlf_lx(nvol), icvrt_lx(nvol), icvlf_ly(nvol), icvrt_ly(nvol), icvlf_lz(nvol), icvrt_lz(nvol))
     allocate(jcvlw_lx(nvol), jcvup_lx(nvol), jcvlw_ly(nvol), jcvup_ly(nvol), jcvlw_lz(nvol), jcvup_lz(nvol))
     allocate(zcvlf_lx(nvol), zcvrt_lx(nvol), zcvlf_ly(nvol), zcvrt_ly(nvol))
     allocate(xld2(nvol), xrd2(nvol), yld2(nvol), yud2(nvol), zld2(nvol), zrd2(nvol))
 
+    write(*,*) 'allocate called'
+
    !  if ((iibm.ne.0).and.(t.ne.0.)) then
    !    xld2(:) = xld(:) + (t-ifirst*dt)*ubcx
    !    xrd2(:) = xrd(:) + (t-ifirst*dt)*ubcx
@@ -79,6 +85,7 @@ contains
     !*****************************************************************
     !! xld,xrd,yld,yud: limits of control volume (!!don't use cex and cey anymore!!)
 
+
     do iv=1,nvol
        ! ok for istret=0 (!!to do for istret=1!!)
        icvlf(iv) = nint(xld(iv)/dx)+1
@@ -288,6 +295,9 @@ contains
     real(mytype), dimension(nz) :: drag1, drag2, drag11, drag22
     real(mytype), dimension(nz) :: drag3, drag4, drag33, drag44
     real(mytype) :: mom1, mom2, mom3, tp1, tp2, tp3, dra1, dra2, dra3
+
+    write(*,*) 'Inside FORCE'
+
   
 
     nvect1=xsize(1)*xsize(2)*xsize(3)
diff --git a/src/ibm.f90 b/src/ibm.f90
index c49a0185..c48ac025 100644
--- a/src/ibm.f90
+++ b/src/ibm.f90
@@ -428,7 +428,7 @@ subroutine cubsplx(u,lind)
   !
   ! Impose the Correct BC
 !   bcimp=lind  
-  write(*,*) lind
+!   write(*,*) lind
   !
   do k=1,xsize(3)
    zm=real(xstart(3)+k-2,mytype)*dz
@@ -457,7 +457,7 @@ subroutine cubsplx(u,lind)
                bcimp=zero
               elseif (lind.eq.1) then
                bcimp=x_pv
-               write(*,*) bcimp
+               ! write(*,*) bcimp
               elseif (lind.eq.2) then
                bcimp=y_pv
               elseif (lind.eq.3) then
diff --git a/src/xcompact3d.f90 b/src/xcompact3d.f90
index f217f993..8585dec7 100644
--- a/src/xcompact3d.f90
+++ b/src/xcompact3d.f90
@@ -239,7 +239,7 @@ subroutine init_xcompact3d()
   endif
 
   if (iforces.eq.1) then
-     call init_forces()
+   !   call init_forces()
      if (irestart==1) then
         call restart_forces(0)
      endif
-- 
GitLab