126
127
logical :: assemble_mass_matrix
127
128
logical :: assemble_inverse_masslump
129
! implicitness parameter, timestep, conservation parameter
130
real :: theta, dt, beta, gravity_magnitude
130
! implicitness parameter, timestep, conservation parameter, nonlinear theta factor
131
real :: theta, dt, beta, gravity_magnitude, itheta
133
! Boundary condition types for velocity, and pressure
134
! the ids have to correspond to the order of the arguments in
135
! the calls to get_entire_boundary_condition below
136
integer, parameter :: BC_TYPE_WEAKDIRICHLET = 1, BC_TYPE_NO_NORMAL_FLOW=2, &
137
BC_TYPE_INTERNAL = 3, BC_TYPE_FREE_SURFACE = 4, &
139
integer, parameter :: PRESSURE_BC_TYPE_WEAKDIRICHLET = 1, PRESSURE_BC_DIRICHLET = 2
132
141
! Stabilisation schemes.
133
142
integer :: stabilisation_scheme
229
238
type(vector_field) :: velocity_bc
230
239
type(scalar_field) :: pressure_bc
231
integer, dimension(:,:), allocatable :: velocity_bc_type, velocity_bc_number
240
integer, dimension(:,:), allocatable :: velocity_bc_type
232
241
integer, dimension(:), allocatable :: pressure_bc_type
234
243
! fields for the assembly of absorption when
489
498
call get_option(trim(u%option_path)//"/prognostic/spatial_discretisation/&
490
499
&conservative_advection", beta)
500
call get_option(trim(u%option_path)//"/prognostic/temporal_discretisation/relaxation", &
492
503
lump_mass=have_option(trim(u%option_path)//&
493
504
&"/prognostic/spatial_discretisation"//&
673
684
ele = fetch(colours(clr), nnid)
674
685
call construct_momentum_element_cg(state, ele, big_m, rhs, ct_m, mass, inverse_masslump, visc_inverse_masslump, &
675
686
x, x_old, x_new, u, oldu, nu, ug, &
677
688
source, absorption, buoyancy, gravity, &
678
689
viscosity, grad_u, &
679
690
tnu, leonard, alpha, &
680
691
gp, surfacetension, &
681
assemble_ct_matrix_here, on_sphere, depth, &
692
assemble_ct_matrix_here, depth, &
682
693
alpha_u_field, abs_wd, temperature, nvfrac, &
683
694
supg_element(thread_num+1))
684
695
end do element_loop
698
709
if((integrate_advection_by_parts.and.(.not.exclude_advection)).or.&
699
710
(integrate_continuity_by_parts)) then
700
711
allocate(velocity_bc_type(u%dim, surface_element_count(u)))
701
allocate(velocity_bc_number(u%dim, surface_element_count(u)))
702
712
call get_entire_boundary_condition(u, &
708
& /), velocity_bc, velocity_bc_type, velocity_bc_number)
719
& /), velocity_bc, velocity_bc_type)
710
721
allocate(pressure_bc_type(surface_element_count(p)))
711
722
call get_entire_boundary_condition(p, &
725
736
! if no_normal flow and no other condition in the tangential directions, or if periodic
726
737
! but not if there's a pressure bc
727
if(((velocity_bc_type(1,sele)==2 .and. sum(velocity_bc_type(:,sele))==2) &
728
.or. any(velocity_bc_type(:,sele)==3)) &
729
.and. pressure_bc_type(sele)==0) cycle
738
if(((velocity_bc_type(1,sele)==BC_TYPE_NO_NORMAL_FLOW &
739
.and. sum(velocity_bc_type(:,sele))==BC_TYPE_NO_NORMAL_FLOW) &
740
.or. any(velocity_bc_type(:,sele)==BC_TYPE_INTERNAL)) &
741
.and. pressure_bc_type(sele)==0) cycle
731
743
ele = face_ele(x, sele)
733
745
call construct_momentum_surface_element_cg(sele, big_m, rhs, ct_m, ct_rhs, &
734
inverse_masslump, x, u, nu, ug, density, p, gravity, &
746
inverse_masslump, x, u, nu, ug, density, gravity, &
735
747
velocity_bc, velocity_bc_type, &
736
748
pressure_bc, pressure_bc_type, &
737
749
assemble_ct_matrix_here, include_pressure_and_continuity_bcs, oldu, nvfrac)
907
919
end subroutine construct_momentum_cg
909
921
subroutine construct_momentum_surface_element_cg(sele, big_m, rhs, ct_m, ct_rhs, &
910
masslump, x, u, nu, ug, density, p, gravity, &
922
masslump, x, u, nu, ug, density, gravity, &
911
923
velocity_bc, velocity_bc_type, &
912
924
pressure_bc, pressure_bc_type, &
913
925
assemble_ct_matrix_here, include_pressure_and_continuity_bcs,&
926
938
type(vector_field), intent(in) :: x, oldu
927
939
type(vector_field), intent(in) :: u, nu
928
940
type(vector_field), pointer :: ug
929
type(scalar_field), intent(in) :: density, p
941
type(scalar_field), intent(in) :: density
930
942
type(vector_field), pointer, intent(in) :: gravity
932
944
type(vector_field), intent(in) :: velocity_bc
944
956
integer :: dim, dim2, i
946
958
integer, dimension(face_loc(u, sele)) :: u_nodes_bdy
947
integer, dimension(face_loc(p, sele)) :: p_nodes_bdy
959
integer, dimension(face_loc(ct_rhs, sele)) :: p_nodes_bdy
948
960
type(element_type), pointer :: u_shape, p_shape
950
962
real, dimension(face_ngi(u, sele)) :: detwei_bdy
951
963
real, dimension(u%dim, face_ngi(u, sele)) :: normal_bdy, upwards_gi
952
real, dimension(u%dim, face_loc(p, sele), face_loc(u, sele)) :: ct_mat_bdy
964
real, dimension(u%dim, face_loc(ct_rhs, sele), face_loc(u, sele)) :: ct_mat_bdy
953
965
real, dimension(u%dim, face_loc(u, sele), face_loc(u, sele)) :: fs_surfacestab
954
966
real, dimension(u%dim, u%dim, face_loc(u, sele), face_loc(u, sele)) :: fs_surfacestab_sphere
955
967
real, dimension(u%dim, u%dim, face_ngi(u, sele)) :: fs_stab_gi_sphere
963
975
real, dimension(u%dim, face_ngi(u, sele)) :: ndotk_k
965
977
u_shape=> face_shape(u, sele)
966
p_shape=> face_shape(p, sele)
978
p_shape=> face_shape(ct_rhs, sele)
968
980
u_nodes_bdy = face_global_nodes(u, sele)
969
p_nodes_bdy = face_global_nodes(p, sele)
981
p_nodes_bdy = face_global_nodes(ct_rhs, sele)
971
983
oldu_val = face_val(oldu, sele)
978
990
! first the advection (dirichlet) bcs:
980
! if no no_normal_flow or free_surface
981
if (velocity_bc_type(1,sele)/=2) then
992
! if no no_normal_flow
993
if (velocity_bc_type(1,sele)/=BC_TYPE_NO_NORMAL_FLOW) then
982
994
if(integrate_advection_by_parts.and.(.not.exclude_advection)) then
984
996
relu_gi = face_val_at_quad(nu, sele)
999
1011
do dim = 1, u%dim
1001
if(velocity_bc_type(dim, sele)==1) then
1013
if(velocity_bc_type(dim, sele)==BC_TYPE_WEAKDIRICHLET) then
1003
1015
call addto(rhs, dim, u_nodes_bdy, -matmul(adv_mat_bdy, &
1004
1016
ele_val(velocity_bc, dim, sele)))
1017
1029
! now do surface integrals for divergence/pressure gradient matrix
1018
1030
if(integrate_continuity_by_parts.and. (assemble_ct_matrix_here .or. include_pressure_and_continuity_bcs)) then
1020
if (velocity_bc_type(1,sele)/=2 .and. velocity_bc_type(1,sele)/=4) then
1032
if (velocity_bc_type(1,sele)/=BC_TYPE_NO_NORMAL_FLOW .and. velocity_bc_type(1,sele)/=BC_TYPE_FREE_SURFACE) then
1022
1034
if(multiphase) then
1023
1035
ct_mat_bdy = shape_shape_vector(p_shape, u_shape, detwei_bdy*face_val_at_quad(nvfrac, sele), normal_bdy)
1030
1042
call addto(ct_rhs, p_nodes_bdy, &
1031
1043
-matmul(ct_mat_bdy(dim,:,:), ele_val(velocity_bc, dim, sele)))
1032
1044
else if (assemble_ct_matrix_here) then
1045
! for open boundaries add in the boundary integral from integrating by parts - for
1046
! other bcs leaving this out enforces a dirichlet-type restriction in the normal direction
1033
1047
call addto(ct_m, 1, dim, p_nodes_bdy, u_nodes_bdy, ct_mat_bdy(dim,:,:))
1035
1049
if(pressure_bc_type(sele)>0) then
1048
1062
! Add free surface stabilisation.
1050
if (velocity_bc_type(1,sele)==4 .and. have_surface_fs_stabilisation) then
1064
if (velocity_bc_type(1,sele)==BC_TYPE_FREE_SURFACE .and. have_surface_fs_stabilisation) then
1051
1065
if (on_sphere) then
1052
1066
upwards_gi=-sphere_inward_normal_at_quad_face(x, sele)
1112
1126
if (assemble_inverse_masslump.and.(.not.(abs_lump_on_submesh))) then
1113
1127
call addto(masslump, u_nodes_bdy, dt*theta*lumped_fs_surfacestab)
1136
if (any(velocity_bc_type(:,sele)==BC_TYPE_FLUX)) then
1138
if(velocity_bc_type(dim,sele)==BC_TYPE_FLUX) then
1139
call addto(rhs, dim, u_nodes_bdy, shape_rhs(u_shape, ele_val_at_quad(velocity_bc, sele, dim)*detwei_bdy))
1122
1145
end subroutine construct_momentum_surface_element_cg
1124
1147
subroutine construct_momentum_element_cg(state, ele, big_m, rhs, ct_m, &
1125
1148
mass, masslump, visc_masslump, &
1126
1149
x, x_old, x_new, u, oldu, nu, ug, &
1128
1151
source, absorption, buoyancy, gravity, &
1129
1152
viscosity, grad_u, &
1130
1153
tnu, leonard, alpha, &
1131
1154
gp, surfacetension, &
1132
assemble_ct_matrix_here, on_sphere, depth, &
1155
assemble_ct_matrix_here, depth, &
1133
1156
alpha_u_field, abs_wd, temperature, nvfrac, supg_shape)
1135
1158
!!< Assembles the local element matrix contributions and places them in big_m
1153
1176
type(vector_field), intent(in) :: x, u, oldu, nu
1154
1177
type(vector_field), pointer :: x_old, x_new, ug
1155
type(scalar_field), intent(in) :: density, p, buoyancy
1178
type(scalar_field), intent(in) :: density, buoyancy, ct_rhs
1156
1179
type(vector_field), intent(in) :: source, absorption, gravity
1157
1180
type(tensor_field), intent(in) :: viscosity, grad_u
1164
1187
type(scalar_field), intent(in) :: gp
1165
1188
type(tensor_field), intent(in) :: surfacetension
1167
logical, intent(in) :: assemble_ct_matrix_here, on_sphere
1190
logical, intent(in) :: assemble_ct_matrix_here
1169
1192
! Wetting and Drying
1170
1193
type(scalar_field), intent(in) :: depth
1190
1213
real, dimension(u%dim, u%dim, ele_ngi(u,ele)) :: J_mat, diff_q
1191
1214
real, dimension(ele_loc(u, ele), ele_ngi(u, ele), u%dim) :: du_t
1192
1215
real, dimension(ele_loc(u, ele), ele_ngi(u, ele), u%dim) :: dug_t
1193
real, dimension(ele_loc(p, ele), ele_ngi(p, ele), u%dim) :: dp_t
1216
real, dimension(ele_loc(ct_rhs, ele), ele_ngi(ct_rhs, ele), u%dim) :: dp_t
1195
1218
real, dimension(u%dim, ele_ngi(u, ele)) :: relu_gi
1196
real, dimension(u%dim, ele_loc(p, ele), ele_loc(u, ele)) :: grad_p_u_mat
1219
real, dimension(u%dim, ele_loc(ct_rhs, ele), ele_loc(u, ele)) :: grad_p_u_mat
1198
1221
! What we will be adding to the matrix and RHS - assemble these as we
1199
1222
! go, so that we only do the calculations we really need
1228
1251
u_ele=>ele_nodes(u, ele)
1229
1252
u_shape=>ele_shape(u, ele)
1231
p_ele=>ele_nodes(p, ele)
1232
p_shape=>ele_shape(p, ele)
1254
p_ele=>ele_nodes(ct_rhs, ele)
1255
p_shape=>ele_shape(ct_rhs, ele)
1234
1257
oldu_val = ele_val(oldu, ele)
1235
1258
! Step 1: Transform
1350
1373
! Buoyancy terms
1351
1374
if(have_gravity) then
1352
call add_buoyancy_element_cg(x, ele, test_function, u, buoyancy, gravity, nvfrac, on_sphere, detwei, rhs_addto)
1375
call add_buoyancy_element_cg(x, ele, test_function, u, buoyancy, gravity, nvfrac, detwei, rhs_addto)
1355
1378
! Surface tension
1683
1706
end subroutine add_sources_element_cg
1685
subroutine add_buoyancy_element_cg(positions, ele, test_function, u, buoyancy, gravity, nvfrac, on_sphere, detwei, rhs_addto)
1708
subroutine add_buoyancy_element_cg(positions, ele, test_function, u, buoyancy, gravity, nvfrac, detwei, rhs_addto)
1686
1709
type(vector_field), intent(in) :: positions
1687
1710
integer, intent(in) :: ele
1688
1711
type(element_type), intent(in) :: test_function
1692
1715
type(scalar_field), intent(in) :: nvfrac
1693
1716
real, dimension(ele_ngi(u, ele)), intent(in) :: detwei
1694
1717
real, dimension(u%dim, ele_loc(u, ele)), intent(inout) :: rhs_addto
1695
logical, intent(in) :: on_sphere
1697
1719
real, dimension(ele_ngi(u, ele)) :: nvfrac_gi
1698
1720
real, dimension(ele_ngi(u, ele)) :: coefficient_detwei
added
removed
assemble/Momentum_CG.F90
