Skip to content

Commit 12b6c56

Browse files
connorairdconnoraird
authored
Merge pull request #17 from UCL-ARC/2-add-src-code
#2 add src code Note: There is an existing bug in the src code. This is fixed in an upcoming PR which add's some test-drive tests
2 parents 2c472fa + 984a053 commit 12b6c56

File tree

9 files changed

+843
-8
lines changed

9 files changed

+843
-8
lines changed

‎.gitignore‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,4 @@
1+
*.o
2+
*.mod
3+
*.DS_Store
4+
build/

‎CMakeLists.txt‎

Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,21 @@
1+
cmake_minimum_required(VERSION 3.9 FATAL_ERROR)
2+
3+
# Set project name
4+
project(
5+
"fortran-tooling"
6+
LANGUAGES "Fortran"
7+
VERSION "0.0.1"
8+
DESCRIPTION "Fortran tooling"
9+
)
10+
11+
set(GFORTRANFLAGS "-ffree-line-length-none")
12+
set(CMAKE_Fortran_FLAGS "${GFORTRAN_FLAGS}") # Update this to match the flags for the chosen compiler
13+
14+
# Define src files
15+
set(PROJ_SRC_DIR "${PROJECT_SOURCE_DIR}/src")
16+
file(GLOB PROJ_MESH_GENERATOR_SOURCES "${PROJ_SRC_DIR}/mesh_generator/*.f90")
17+
file(GLOB PROJ_POISSON_SOURCES "${PROJ_SRC_DIR}/poisson/*.f90")
18+
19+
# Build src executables
20+
add_executable("${PROJECT_NAME}-mesh-generator" "${PROJ_MESH_GENERATOR_SOURCES}")
21+
add_executable("${PROJECT_NAME}-poisson" "${PROJ_POISSON_SOURCES}")

‎README.md‎

Lines changed: 33 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -25,3 +25,36 @@ This repository aims to improve Fortran best practices within UCL and the wider
2525
| FruitDemand | Apr 2021 | Mar 2023 | | [#382](https://github.com/UCL-ARC/research-software-opportunities/issues/382) | Make, Ford, PFUnit | |
2626
| Trove | Jan 2021 | Aug 2021 | | [#404](https://github.com/UCL-ARC/research-software-opportunities/issues/404) | | |
2727
| Zacros | Jan 2021 | Sep 2022 | | [#349](https://github.com/UCL-ARC/research-software-opportunities/issues/349) & older | CMake, CTest | |
28+
29+
## src code
30+
31+
There are two src codes within this repository [mesh_generator](./src/mesh_generator/) and [poisson](./src/poisson/). These are designed to work together.
32+
- `mesh_generator` generates a basic square 2D triangular mesh (see [mesh_generator.f90](./src/mesh_generator/mesh_generator.f90) for more details).
33+
- `poisson` is a solver which finds the solution of the steady-state heat conduction equation represented by the Poisson equation over a 2D traingular mesh (see [poisson.f90](./src/poisson/poisson.f90) for more details).
34+
35+
## Building
36+
37+
We are utilising cmake (see [CMakeLists.txt](./CMakeLists.txt)) for our build system. Therefore, to build this repository, please run the following
38+
```sh
39+
cmake -B build
40+
```
41+
This will create a [build](./build) directory from within which the project can be compiled...
42+
```sh
43+
cd build
44+
make
45+
```
46+
This will produce executables for the two src codes, `fortran-tooling-mesh-generator` and `fortran-tooling-poisson`.
47+
48+
## Running the src
49+
50+
### Mesh generator
51+
52+
```sh
53+
./build/fortran-tooling-mesh-generator <box_size> <edge_size>
54+
```
55+
56+
### Poisson solver
57+
58+
```sh
59+
./build/fortran-tooling-poisson # then respond to prompt with the mesh name, likely to be `square_mesh`
60+
```

‎src/mesh_generator/main.f90‎

Lines changed: 44 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,44 @@
1+
program main
2+
use, intrinsic :: iso_fortran_env
3+
use mesh_generator
4+
5+
implicit none
6+
7+
! Command line arguments
8+
integer :: argl
9+
character(len=:), allocatable :: a
10+
integer(kind=int64) :: box_size
11+
real(kind=real64) :: edge_size
12+
13+
! Mesh variables
14+
integer(kind=int64) :: num_nodes, num_elements, num_boundary_nodes, &
15+
num_edges_per_boundary
16+
integer(kind=int64), dimension(:, :), allocatable :: elements, boundary_edges
17+
real(kind=real64), dimension(:, :), allocatable :: nodes
18+
19+
! Get command line args (Fortran 2003 standard)
20+
if (command_argument_count() > 0) then
21+
call get_command_argument(1, length=argl)
22+
allocate(character(argl) :: a)
23+
call get_command_argument(1, a)
24+
read(a,*) box_size
25+
call get_command_argument(2, a)
26+
read(a,*) edge_size
27+
end if
28+
29+
! Output start message
30+
write(*,'(A)') "Generating mesh using:"
31+
write(*,'(A,1I16)') "box size: ", box_size
32+
write(*,*) " process: ", edge_size
33+
34+
call calculate_mesh_parameters(box_size, edge_size, num_edges_per_boundary, num_nodes, num_boundary_nodes, num_elements)
35+
36+
! Allocate arrays
37+
allocate(nodes(2, num_nodes))
38+
allocate(elements(3, num_elements))
39+
allocate(boundary_edges(3, num_boundary_nodes))
40+
41+
call calculate_mesh(num_edges_per_boundary, num_nodes, num_elements, num_boundary_nodes, nodes, elements, boundary_edges)
42+
43+
call write_mesh_to_file(num_nodes, num_elements, num_boundary_nodes, nodes, elements, boundary_edges)
44+
end program main

‎src/mesh_generator/mesh_generator.f90‎

Lines changed: 185 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,185 @@
1+
!!----------------------------------------------------------------------------*
2+
!! *** Program MESH_GENERATOR *** *
3+
!! *
4+
!! This module generates a basic square 2D traingular mesh. The *
5+
!! geometry of this mesh is kept simple by placing nodes in a *
6+
!! regular pattern, as demonstrated below (note the arrangment *
7+
!! of the node and element IDs). *
8+
!! *
9+
!! 21._____22._____23._____24._____25. *
10+
!! | / | / | / | / | *
11+
!! | 29/ 25| 30/ 26| 31/ 27| 32/ 28| *
12+
!! | / | / | / | / | *
13+
!! 16._____17._____18._____19._____20. *
14+
!! | / | / | / | / | *
15+
!! | 21/ 17| 22/ 18| 23/ 19| 24/ 20| *
16+
!! | / | / | / | / | *
17+
!! 11._____12._____13._____14._____15. *
18+
!! | / | / | / | / | *
19+
!! | 13/ 9| 14/ 10| 15/ 11| 16/ 12| *
20+
!! | / | / | / | / | *
21+
!! 6.______7.______8.______9._____10. *
22+
!! | / | / | / | / | *
23+
!! | 5/ 1| 6/ 2| 7/ 3| 8/ 4| *
24+
!! | / | / | / | / | *
25+
!! 1.______2.______3.______4.______5. *
26+
!! *
27+
!! The size of the output square and its elements is determined *
28+
!! by two inputs, box_size (the length of each side of the *
29+
!! square) and edge_size (the length of a single elements *
30+
!! boundary edge). Therefore, the above mesh would be generated *
31+
!! by inputs that satisfy - box_size / edge_size = 4. *
32+
!! *
33+
!!----------------------------------------------------------------------------*
34+
module mesh_generator
35+
use, intrinsic :: iso_fortran_env
36+
implicit none
37+
38+
contains
39+
40+
subroutine calculate_mesh_parameters(box_size, edge_size, num_edges_per_boundary, num_nodes, num_boundary_nodes, num_elements)
41+
implicit none
42+
integer(kind=int64), intent(in) :: box_size
43+
real(kind=real64), intent(in) :: edge_size
44+
integer(kind=int64), intent(out) :: num_edges_per_boundary, num_nodes, num_boundary_nodes, num_elements
45+
46+
num_edges_per_boundary = floor(box_size / edge_size)
47+
num_nodes = (num_edges_per_boundary + 1)**2
48+
num_boundary_nodes = (num_edges_per_boundary) * 4
49+
num_elements = 2 * (num_edges_per_boundary)**2
50+
51+
! write(*,*) "** Mesh Parameters **"
52+
! write(*,*) " num_edges_per_boundary:", num_edges_per_boundary
53+
! write(*,*) " num_nodes:", num_nodes
54+
! write(*,*) " num_boundary_nodes:", num_boundary_nodes
55+
! write(*,*) " num_elements:", num_elements
56+
end subroutine calculate_mesh_parameters
57+
58+
subroutine calculate_mesh(num_edges_per_boundary, num_nodes, num_elements, num_boundary_nodes, nodes, elements, boundary_edges)
59+
implicit none
60+
integer(kind=int64), intent(in) :: num_edges_per_boundary, num_nodes, num_boundary_nodes, num_elements
61+
integer(kind=int64), dimension(3, num_elements), intent(inout) :: elements
62+
integer(kind=int64), dimension(3, num_boundary_nodes), intent(inout) :: boundary_edges
63+
real(kind=real64), dimension(2, num_nodes), intent(inout) :: nodes
64+
65+
integer :: num_nodes_per_boundary, bottom_left_node, counter, i, j
66+
67+
num_nodes_per_boundary = num_edges_per_boundary + 1
68+
69+
counter = 1
70+
do i = 1, num_nodes_per_boundary
71+
do j = 1, num_nodes_per_boundary
72+
nodes(1, counter) = i
73+
nodes(2, counter) = j
74+
counter = counter + 1
75+
end do
76+
end do
77+
78+
counter = 1
79+
do i = 1, num_edges_per_boundary
80+
do j = 1, num_edges_per_boundary
81+
bottom_left_node = j + (i - 1) * num_nodes_per_boundary
82+
83+
elements(1, counter) = bottom_left_node
84+
elements(2, counter) = bottom_left_node + 1
85+
elements(3, counter) = bottom_left_node + 1 + num_nodes_per_boundary
86+
87+
elements(1, counter + num_edges_per_boundary) = bottom_left_node
88+
elements(2, counter + num_edges_per_boundary) = bottom_left_node + num_nodes_per_boundary + 1
89+
elements(3, counter + num_edges_per_boundary) = bottom_left_node + num_nodes_per_boundary
90+
91+
counter = counter + 1
92+
end do
93+
counter = counter + num_edges_per_boundary
94+
end do
95+
96+
! If we are along the bottom boundary
97+
do i = 1, num_edges_per_boundary
98+
! bottom boundary
99+
boundary_edges(1, i) = i ! left node
100+
boundary_edges(2, i) = i + 1 ! right node
101+
boundary_edges(3, i) = i*2 - 1 ! element
102+
103+
! right boundary
104+
boundary_edges(1, i + num_edges_per_boundary) = i * num_nodes_per_boundary
105+
boundary_edges(2, i + num_edges_per_boundary) = (i + 1) * num_nodes_per_boundary
106+
boundary_edges(3, i + num_edges_per_boundary) = (2*i - 1) * num_edges_per_boundary
107+
108+
! top boundary
109+
boundary_edges(1, i + num_edges_per_boundary * 2) = num_nodes - i + 1
110+
boundary_edges(2, i + num_edges_per_boundary * 2) = num_nodes - i
111+
boundary_edges(2, i + num_edges_per_boundary * 2) = num_elements - i + 1
112+
113+
! left boundary
114+
boundary_edges(1, i + num_edges_per_boundary * 3) = (num_nodes_per_boundary - i) * num_nodes_per_boundary + 1
115+
boundary_edges(2, i + num_edges_per_boundary * 3) = (num_nodes_per_boundary - 1 - i) * num_nodes_per_boundary + 1
116+
boundary_edges(3, i + num_edges_per_boundary * 3) = num_elements - (num_edges_per_boundary - 1) - (2 * (i - 1) * num_edges_per_boundary)
117+
end do
118+
119+
end subroutine calculate_mesh
120+
121+
subroutine write_mesh_to_file(num_nodes, num_elements, num_boundary_nodes, nodes, elements, boundary_edges)
122+
implicit none
123+
integer(kind=int64), intent(in) :: num_nodes, num_elements, num_boundary_nodes
124+
integer(kind=int64), dimension(3, num_elements), intent(inout) :: elements
125+
integer(kind=int64), dimension(3, num_boundary_nodes), intent(inout) :: boundary_edges
126+
real(kind=real64), dimension(2, num_nodes), intent(inout) :: nodes
127+
128+
character*11 :: file_name
129+
integer :: file_io
130+
integer :: iostat
131+
integer :: i, num_sets, num_dirichlet_boundary_conditions, num_neumann_boundary_conditions
132+
133+
! Baked in defaults
134+
num_sets = 1
135+
num_dirichlet_boundary_conditions = 1
136+
num_neumann_boundary_conditions = 0
137+
138+
file_name = "square_mesh"
139+
file_io = 100
140+
141+
! Write outpout
142+
open (unit=file_io, &
143+
file=file_name, &
144+
status="replace", &
145+
IOSTAT=iostat)
146+
147+
if( iostat .ne. 0) then
148+
write(*,'(a)') ' *** Error when opening '//trim(file_name)
149+
else
150+
write(*,'(/,a)') ' *** '//trim(file_name)//' opened'
151+
end if
152+
153+
write(file_io,*) "! num_nodes, num_elements, num_boundary_points, num_sets, num_dirichlet_boundary_conditions, num_neumann_boundary_conditions"
154+
write(file_io,*) num_nodes, num_elements, num_boundary_nodes, num_sets, num_dirichlet_boundary_conditions, num_neumann_boundary_conditions
155+
156+
write(file_io,*) "! jb,vb(1,jb),vb(2,jb),vb(3,jb) - as many lines as num_sets"
157+
write(file_io,*) 1, 1, 1, 1
158+
159+
write(file_io,*) "! jb,vb1(jb) - as many lines as num_dirichlet_boundary_conditions"
160+
write(file_io,*) 1, 0
161+
162+
write(file_io,*) "! jb,vb2(jb) - as many lines as num_neumann_boundary_conditions"
163+
164+
write(file_io,*) "! jp,coordinates(1,jp),coordinates(2,jp) - as many lines as num_nodes"
165+
do i = 1, num_nodes
166+
write(file_io,*) i, nodes(1, i), nodes(2, i)
167+
end do
168+
169+
write(file_io,*) "! je,element_to_node(1,je),element_to_node(2,je),element_to_node(3,je),vb_index(je) - as many lines as num_elements"
170+
do i = 1, num_elements
171+
write(file_io,*) i, elements(1, i), elements(2, i), elements(3, i), 1
172+
end do
173+
174+
write(file_io,*) "! boundary_node_num(1,ib),boundary_node_num(2,ib) - as many lines as num_boundary_points"
175+
do i = 1, num_boundary_nodes
176+
write(file_io,*) i, 1
177+
end do
178+
179+
write(file_io,*) "! num_side_nodes(1,ib),num_side_nodes(2,ib),num_side_nodes(3,ib),num_side_nodes(4,ib) - as many lines as num_boundary_points"
180+
do i = 1, num_boundary_nodes
181+
write(file_io,*) boundary_edges(1, i), boundary_edges(2, i), boundary_edges(3, i), 0
182+
end do
183+
end subroutine write_mesh_to_file
184+
185+
end module mesh_generator

‎src/poisson/main.f90‎

Lines changed: 44 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,44 @@
1+
program main
2+
use poisson
3+
4+
implicit none
5+
6+
integer :: boundary_index(mxp), element_to_node(3,mxe), &
7+
vb_index(mxe), boundary_node_num(2,mxb), num_side_nodes(4,mxb)
8+
real :: coordinates(2, mxp), nodal_value_of_f(mxp), rhs_vector(mxp), beta(mxp), f_increment(mxp), &
9+
vb1(mxc), vb2(mxc), vb(3,mxc), element_stiffness(6,mxe), &
10+
pre_conditioning_matrix(mxp)
11+
integer :: fname_io = 100, fname_out_io = 101
12+
character*120 :: input_fname, output_fname
13+
14+
!!
15+
!! *** Asks for I/O file names
16+
!!
17+
write(*,'(10(/),7(a,/),//)') &
18+
' *********************************', &
19+
' *** ***', &
20+
' *** * P O I S S O N * ***', &
21+
' *** ***', &
22+
' *** PCG-FEM HEAT CONDUCTION ***', &
23+
' *** ***', &
24+
' *********************************'
25+
input_fname = textread(' Enter input file name : ')
26+
output_fname = trim(input_fname)//'.out'
27+
call open_file(input_fname, 'old', fname_io)
28+
call open_file(output_fname, 'new', fname_out_io)
29+
30+
!!
31+
!! *** Reads the triangular mesh and problem constants: Kx,Ky,Q,fp,q
32+
!!
33+
call inp(element_to_node,vb_index,coordinates,boundary_node_num,num_side_nodes,vb,vb1,vb2,fname_io)
34+
35+
!!
36+
!! *** Assembles and solves the system of equations
37+
!!
38+
call pcg(element_to_node,vb_index,coordinates,nodal_value_of_f,boundary_node_num,num_side_nodes,vb,vb1,vb2,element_stiffness,rhs_vector,beta,f_increment,boundary_index,pre_conditioning_matrix)
39+
40+
!!
41+
!! *** Writes the computed solution
42+
!!
43+
call out(element_to_node,coordinates,nodal_value_of_f,fname_out_io)
44+
end program main

0 commit comments

Comments
 (0)