Add files via upload

Documents/steps_recap.txt

@@ -0,0 +1,7 @@
+In u_star and v_star calculation
+
+Velocities with index 1/2 are represented by letters (8 in total). This is due to the conditions on u_tilde et v_tilde.
+
+(a,b,c,d) for calculations in x
+(e,f,g,h) for calculations in y
+

Makefile

@@ -0,0 +1,78 @@
+# Makefile 
+# Nom du compilateur
+FC = gfortran
+
+# Options de compilation: optimisation, debug etc...
+OPT = -fdefault-real-8 -O3
+# nom de l'executable
+EXE = NS_lid_driven_cavity
+# Options de l'edition de lien..
+LINKOPT = 
+
+# Defining the objects (OBJS) variables
+
+OBJS =  \
+       variables_module.o \
+       main.o \
+       reading_data.o \
+       mesh_creation.o \
+       initial_conditions.o \
+       boundary_conditions.o \
+       timestep.o \
+       discretisation_upwind.o \
+       discretisation_centre.o \
+       calcul_2nd_membre.o \
+       matgen_cavity.o \
+       solver.o \
+       calcul_pressure_field.o \
+       calcul_new_velocity.o \
+       velocities_output.o \
+       ensight_cavity.o \
+
+
+
+# Linking object files
+exe :   $(OBJS)
+	$(FC) $(LINKOPT) $(OBJS)  -o $(EXE) 
+
+%.o:%.f90
+	$(FC) $(OPT) -c $<
+
+main.o : variables_module.o
+
+reading_data.o : variables_module.o
+
+mesh_creation.o: variables_module.o
+
+initial_conditions.o : variables_module.o
+
+boundary_conditions.o : variables_module.o
+
+timestep.o : variables_module.o
+
+discretisation_upwind.o : variables_module.o
+
+discretisation_centre.o : variables_module.o
+
+calcul_2nd_membre.o : variables_module.o
+
+calcul_pressure_field.o : variables_module.o
+
+calcul_new_velocity.o : variables_module.o
+
+velocities_output.o : variables_module.o
+
+
+
+# Removing object files
+clean :
+	/bin/rm -f $(OBJS) $(EXE)  *.mod
+
+cleanall : 
+	/bin/rm -f $(OBJS) $(EXE)  *.mod
+	/bin/rm -f *.scl
+	/bin/rm -f *.vec
+
+config :
+	if [ ! -d obj ] ; then mkdir obj ; fi
+	if [ ! -d run ] ; then mkdir run ; fi

boundary_conditions.f90

@@ -0,0 +1,77 @@
+subroutine boundary_conditions()
+use variables
+implicit none
+
+
+!---------------------------------------------------!
+!         Boundary conditions calculation           !
+!---------------------------------------------------!
+
+
+!!!!!!!!Velocity boundary conditions!!!!!!!!
+
+!Left vertical wall
+do j=0,ny+1
+  u(0,j) = 0
+  v(0,j) = -v(1,j)
+  u1(0,j) = 0
+  v1(0,j) = -v1(1,j)
+  u_star(0,j) = 0 
+  v_star(0,j) = -v_star(1,j)
+end do
+
+!Right vertical wall
+do j=0,ny+1
+  u(nx,j) = 0
+  v(nx+1,j) = -v(nx,j)
+  u1(nx,j) = 0
+  v1(nx+1,j) = -v1(nx,j)
+  u_star(nx,j) = 0
+  v_star(nx+1,j) = -v_star(nx,j)
+end do
+
+!Bottom horizontal wall
+do i=0,nx+1
+  u(i,0) = -u(i,1)
+  v(i,0) = 0
+  u1(i,0) = -u1(i,1)
+  v1(i,0) = 0
+  u_star(i,0) = -u_star(i,1)
+  v_star(i,0) = 0
+end do
+
+!Top horizontal wall
+do i=0,nx+1
+  u(i,ny+1) = 2-u(i,ny) 
+  v(i,ny) = 0
+  u1(i,ny+1) = 2-u1(i,ny)
+  v1(i,ny) = 0
+  u_star(i,ny+1) = 2-u_star(i,ny)
+  v_star(i,ny) = 0
+end do
+
+
+!!!!!!!Pressure boundary conditions!!!!!!!
+
+!Left vertical wall
+do j=0,ny+1
+   p(0,j)=p(1,j)
+end do
+
+!Right vertical wall
+do j=0,ny+1
+   p(nx+1,j)=p(nx,j)
+end do
+
+
+!Bottom horizontal wall
+do i=0,nx+1
+   p(i,0)=p(i,1)
+end do
+
+!Top horizontal wall
+do i=0,nx+1
+   p(i,ny+1)=p(i,ny)
+end do
+
+end subroutine boundary_conditions

calcul_2nd_membre.f90

@@ -0,0 +1,30 @@
+subroutine calcul_2nd_membre() ! Calculate the second membre of equation AX=B
+use variables
+implicit none
+
+
+	do j=1,ny
+	   do i=1,nx
+	      div(i,j)= ddt*(ddx*(u_star(i,j)-u_star(i-1,j))+ddy*(v_star(i,j)-v_star(i,j-1))) 
+	   enddo
+	enddo
+
+	sum=0.0
+	do j=1,ny
+	   do i=1,nx
+	      sum=sum+div(i,j)
+	   enddo
+	enddo
+
+        sum=sum/float(nx*ny)
+	k=1
+	do j=1,ny
+	   do i=1,nx
+	      div(i,j)=div(i,j)-sum
+              div1(k)=-div(i,j)
+	      x1(k)=0.
+	      k=k+1
+	    enddo
+	enddo
+
+end subroutine calcul_2nd_membre

calcul_new_velocity.f90

@@ -0,0 +1,87 @@
+subroutine calcul_new_velocity()
+use variables
+implicit none
+
+
+!-------------------------------------------------!
+!    Calculation of velocity field at t = dt*n+1  !
+!-------------------------------------------------!
+
+!!! In x direction !!!
+
+do j=1,ny
+  do i=1,nx-1
+u1(i,j) = u_star(i,j) - dt*ddx*(p(i+1,j)-p(i,j))
+  end do
+end do
+
+!!! In y direction !!!
+
+do j=1,ny-1
+  do i =1,nx
+v1(i,j) = v_star(i,j) - dt*ddy*(p(i,j+1)-p(i,j))
+  end do
+end do
+
+
+
+
+!------------------------------------------------------------------------------------------!
+!               Centering variables to get matrix (1:nx,1:ny) for Paraview                 !
+!------------------------------------------------------------------------------------------!
+
+
+do j=1,ny
+  do i=1,nx
+uc(i,j) = 0.5*(u1(i,j)+u1(i-1,j))
+vc(i,j) = 0.5*(v1(i,j)+v1(i,j-1))
+  end do
+end do
+
+
+
+!---------------------------------------------------!
+!    Rotational calculus (corners + centering)      !
+!---------------------------------------------------!
+
+
+!!! Corner boundary conditions (variable w) !!!
+
+w(1,1)=0
+w(nx,1)=0
+w(1,ny)=0
+w(nx,ny)=0
+
+!!! General calculation !!!
+
+do j=2,ny
+  do i=2,nx
+w(i,j) = ddx*(v1(i+1,j)-v1(i,j))-ddy*(u1(i,j+1)-u1(i,j))
+  end do
+end do
+
+!!! Top and bottom wall computation !!!
+
+do i=2,nx-1
+  w(i,1) = ddx*(v1(i+1,1)-v1(i,1))-ddy*(u1(i,2)-u1(i,1))
+  w(i,ny) = ddx*(v1(i+1,ny)-v1(i,ny))-ddy*(u1(i,ny+1)-u1(i,j))
+end do
+
+!!! Left and right wall computation !!!
+
+do j=2,ny-1
+  w(1,j) = ddx*(v1(2,j)-v1(1,j))-ddy*(u1(1,j+1)-u1(1,j))
+  w(nx,j) = ddx*(v1(nx+1,j)-v1(nx,j))-ddy*(u1(nx,j+1)-u1(nx,j))
+end do
+
+!!! Centered rotational (at the center of mesh cells) !!!
+
+do j=2,ny
+  do i=2,nx
+rot(i,j) = 0.25*(w(i,j)+w(i-1,j)+w(i,j-1)+w(i-1,j-1))
+  end do
+end do
+
+end subroutine calcul_new_velocity
+
+

calcul_pressure_field.f90

@@ -0,0 +1,25 @@
+subroutine calcul_pressure_field()
+use variables
+implicit none
+
+!-------------------------------------------------!
+!              Update pressure field              !
+!-------------------------------------------------!
+
+    k=1
+       do j=1,ny 
+	  do i=1,nx 
+  	      pre(i,j)=x1(k)
+	      k=k+1
+	  enddo
+       enddo
+
+
+       do j=1,ny
+          do i=1,nx
+             p(i,j) = pre(i,j)
+          end do
+       end do
+
+
+end subroutine calcul_pressure_field