« back to all changes in this revision
Viewing changes to assemble/Hybridized_Helmholtz.F90
- browse files at revision 3602
- compare with another revision
- download diff
- download tarball
- view history from revision 3602
- Committer: Jemma Shipton
- Date: 2012-10-23 12:39:04 UTC
- mfrom: (3574.2.151 bdfm1-nonlinear-sw)
- Revision ID: jshipton@ic.ac.uk-20121023123904-oalpymh8vx716nj6
merging changes from Colin's bdfm1-nonlinear-sw branch
- files added:
- assemble/Bubble_tools.F90
- tests/bdfm1-nonlinear
- tests/galewsky
- tests/lerouxeddy
- tests/sphere_williamson2
- tests/sphere_williamson3
- tests/sphere_williamson5
- tests/sw_advection_tests
- tests/sw_sphere_advection
- files modified:
- Makefile.in
added
removed
assemble/Hybridized_Helmholtz.F90
24
24
! License along with this library; if not, write to the Free Software
25
25
! Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
26
26
! USA
27
28
27
29
28
#include "fdebug.h"
30
29
49
48
use manifold_tools
50
49
implicit none
51
50
51
!Variables belonging to the Newton solver
52
logical, private :: newton_initialised = .false.
53
!Cached Newton solver matrices
54
type(csr_matrix), private :: Newton_lambda_mat
55
type(block_csr_matrix), private :: Newton_continuity_mat
56
type(real_matrix), pointer, dimension(:), private &
57
&:: Newton_local_solver_cache=>null()
58
type(real_matrix), pointer, dimension(:), private &
59
&:: Newton_local_solver_rhs_cache=>null()
60
52
61
contains
62
63
subroutine solve_hybridised_timestep_residual(state,newU,newD,&
64
UResidual, DResidual)
65
66
! Subroutine to apply one Newton iteration to hybridised shallow
67
! water equations
68
! Written to cache all required matrices
69
implicit none
70
type(state_type), intent(inout) :: state
71
type(vector_field), intent(inout) :: newU !U at next timestep
72
type(scalar_field), intent(inout) :: newD !D at next timestep
73
type(vector_field), intent(in), optional :: UResidual
74
type(scalar_field), intent(in), optional :: DResidual
75
!
76
type(vector_field), pointer :: X, U, down, U_cart, U_old
77
type(scalar_field), pointer :: D,f, D_old
78
type(scalar_field) :: lambda, D_res, D_res2
79
type(vector_field) :: U_res, U_res2
80
type(scalar_field), target :: lambda_rhs, u_cpt
81
type(csr_sparsity) :: lambda_sparsity, continuity_sparsity
82
type(csr_matrix) :: continuity_block_mat
83
type(mesh_type), pointer :: lambda_mesh
84
real :: D0, dt, g, theta, tolerance
85
integer :: ele,i1, stat, dim1, dloc, uloc,mdim,n_constraints
86
logical :: have_constraint
87
character(len=OPTION_PATH_LEN) :: constraint_option_string
88
89
ewrite(1,*) 'solve_hybridised_timestep_residual(state)'
90
91
ewrite(1,*) 'TIMESTEPPING LAMBDA'
92
93
!get parameters
94
call get_option("/physical_parameters/gravity/magnitude", g)
95
call get_option("/timestepping/theta",theta)
96
call get_option("/material_phase::Fluid/scalar_field::LayerThickness/&
97
&prognostic/mean_layer_thickness",D0)
98
call get_option("/timestepping/timestep", dt)
99
100
!Pull the fields out of state
101
D=>extract_scalar_field(state, "LayerThickness")
102
f=>extract_scalar_field(state, "Coriolis")
103
U=>extract_vector_field(state, "LocalVelocity")
104
X=>extract_vector_field(state, "Coordinate")
105
down=>extract_vector_field(state, "GravityDirection")
106
U_cart => extract_vector_field(state, "Velocity")
107
U_old => extract_vector_field(state, "OldLocalVelocity")
108
D_old => extract_scalar_field(state, "OldLayerThickness")
109
110
!Allocate local field variables
111
lambda_mesh=>extract_mesh(state, "VelocityMeshTrace")
112
call allocate(lambda,lambda_mesh,name="LagrangeMultiplier")
113
call allocate(lambda_rhs,lambda%mesh,"LambdaRHS")
114
call zero(lambda_rhs)
115
call allocate(U_res,U%dim,U%mesh,"VelocityResidual")
116
call allocate(D_res,D%mesh,"LayerThicknessResidual")
117
call allocate(U_res2,U%dim,U%mesh,"VelocityResidual2")
118
call allocate(D_res2,D%mesh,"LayerThicknessResidual2")
119
120
if(.not.newton_initialised) then
121
!construct/extract sparsities
122
lambda_sparsity=get_csr_sparsity_firstorder(&
123
&state, lambda%mesh, lambda&
124
&%mesh)
125
continuity_sparsity=get_csr_sparsity_firstorder(state, u%mesh,&
126
& lambda%mesh)
127
128
!allocate matrices
129
call allocate(Newton_lambda_mat,lambda_sparsity)
130
call zero(Newton_lambda_mat)
131
call allocate(Newton_continuity_mat,continuity_sparsity,(/U%dim,1/))
132
call zero(Newton_continuity_mat)
133
134
mdim = mesh_dim(U)
135
have_constraint = &
136
&have_option(trim(U%mesh%option_path)//"/from_mesh/constraint_type")
137
138
allocate(newton_local_solver_cache(ele_count(U)))
139
allocate(newton_local_solver_rhs_cache(ele_count(U)))
140
do ele = 1, ele_count(D)
141
uloc = ele_loc(U,ele)
142
dloc = ele_loc(d,ele)
143
n_constraints = 0
144
if(have_constraint) then
145
n_constraints = ele_n_constraints(U,ele)
146
allocate(newton_local_solver_cache(ele)%ptr(&
147
mdim*uloc+dloc+n_constraints,&
148
mdim*uloc+dloc+n_constraints))
149
allocate(newton_local_solver_rhs_cache(ele)%ptr(&
150
mdim*uloc+dloc,mdim*uloc+dloc))
151
end if
152
end do
153
154
!Assemble matrices
155
do ele = 1, ele_count(D)
156
call assemble_newton_solver_ele(D,f,U,X,down,lambda_rhs,ele, &
157
&g,dt,theta,D0,&
158
&lambda_mat=Newton_lambda_mat,&
159
&continuity_mat=Newton_continuity_mat,&
160
&Local_solver_matrix=Newton_local_solver_cache(ele)%ptr,&
161
&Local_solver_rhs=Newton_local_solver_rhs_cache(ele)%ptr)
162
end do
163
newton_initialised = .true.
164
end if
165
166
if(present(DResidual).and.present(UResidual).and..not.&
167
have_option('/material_phase::Fluid/vector_field::Velocity/prognostic/wave_equation/fully_coupled/linear_debug')) then
168
!Set residuals from nonlinear input
169
call set(U_res,UResidual)
170
call set(D_res,DResidual)
171
else
172
!Compute residuals for linear equation
173
do ele = 1, ele_count(D)
174
call get_linear_residuals_ele(U_res,D_res,&
175
U_old,D_old,newU,newD,&
176
newton_local_solver_cache(ele)%ptr,&
177
newton_local_solver_rhs_cache(ele)%ptr,ele)
178
end do
179
ewrite(2,*) 'cjc U_res calc', sum(newU%val), maxval(abs(U_res%val))
180
ewrite(2,*) 'cjc D_res calc', sum(newD%val), maxval(abs(D_res%val))
181
end if
182
183
do ele = 1, ele_count(D)
184
call local_solve_residuals_ele(U_res,D_res,&
185
newton_local_solver_cache(ele)%ptr,ele)
186
end do
187
188
call mult_t(lambda_rhs,Newton_continuity_mat,U_res)
189
call scale(lambda_rhs,-1.0)
190
ewrite(2,*) 'LAMBDARHS', maxval(abs(lambda_rhs%val))
191
192
call zero(lambda)
193
!Solve the equations
194
call petsc_solve(lambda,newton_lambda_mat,lambda_rhs,&
195
option_path=trim(U_cart%mesh%option_path)//&
196
&"/from_mesh/constraint_type")
197
ewrite(2,*) 'LAMBDA', maxval(abs(lambda%val))
198
199
!Update new U and new D from lambda
200
!Compute residuals
201
if(present(DResidual).and.present(UResidual).and..not.&
202
have_option('/material_phase::Fluid/vector_field::Velocity/prognostic/wave_equation/fully_coupled/linear_debug')) then
203
call set(U_res,UResidual)
204
call set(D_res,DResidual)
205
else
206
do ele = 1, ele_count(D)
207
call get_linear_residuals_ele(U_res,D_res,&
208
U_old,D_old,newU,newD,&
209
newton_local_solver_cache(ele)%ptr,&
210
newton_local_solver_rhs_cache(ele)%ptr,ele)
211
end do
212
end if
213
ewrite(2,*) 'cjc U_res recalc', sum(newU%val), maxval(abs(U_res%val))
214
ewrite(2,*) 'cjc D_res recalc', sum(newD%val), maxval(abs(D_res%val))
215
216
call mult(U_res2,Newton_continuity_mat,lambda)
217
call addto(U_res,U_res2)
218
do ele = 1, ele_count(U)
219
call local_solve_residuals_ele(U_res,D_res,&
220
newton_local_solver_cache(ele)%ptr,ele)
221
end do
222
223
!Negative scaling occurs here.
224
call scale(U_res,-1.0)
225
call scale(D_res,-1.0)
226
227
ewrite(2,*) 'Delta U', maxval(abs(U_res%val))
228
ewrite(2,*) 'Delta D', maxval(abs(D_res%val))
229
230
call addto(newU,U_res)
231
call addto(newD,D_res)
232
233
!Deallocate local variables
234
call deallocate(lambda_rhs)
235
call deallocate(lambda)
236
call deallocate(U_res)
237
call deallocate(D_res)
238
call deallocate(U_res2)
239
call deallocate(D_res2)
240
241
ewrite(1,*) 'END solve_hybridised_timestep_residual(state)'
242
243
end subroutine solve_hybridised_timestep_residual
244
53
245
subroutine solve_hybridized_helmholtz(state,D_rhs,U_Rhs,&
54
246
&D_out,U_out,&
55
247
&dt_in,theta_in,&
65
257
! <<\gamma, [u]>> = 0
66
258
! (i.e. for unit testing)
67
259
! otherwise solve:
68
! <w,u> + dt*theta*<w,fu^\perp> - dt*theta*g <div w,d> + <<[w],d>> =
260
! <w,u> + dt*theta*<w,fu^\perp> - dt*theta*g <div w,d> + <<[w],l>> =
69
261
! -dt*<w f(u^n)^\perp> + dt*g*<div w, d^n>
70
262
! <\phi,\eta> + dt*theta*<\phi,div u> = <\ph
71
263
! <<\gamma, [u]>> = 0
253
445
end if
254
446
end if
255
447
256
if(have_option('/geometry/mesh::VelocityMesh/from_mesh/constraint_type/c&
257
&heck_continuity')) then
448
if(have_option('/geometry/mesh::VelocityMesh/check_continuity')) then
258
449
U_cart => extract_vector_field(state, "Velocity")
259
450
if(present(U_out)) then
260
451
call project_local_to_cartesian(X,U_out,U_cart,weights=weights)
268
459
end if
269
460
270
461
call deallocate(lambda_mat)
462
call deallocate(continuity_mat)
271
463
call deallocate(lambda_rhs)
272
464
call deallocate(lambda)
273
465
deallocate(weights)
275
467
ewrite(1,*) 'END subroutine solve_hybridized_helmholtz'
276
468
277
469
end subroutine solve_hybridized_helmholtz
470
471
subroutine get_linear_residuals_ele(U_res,D_res,U,D,&
472
newU,newD,local_solver,local_solver_rhs,ele)
473
!Subroutine
474
type(vector_field), intent(inout) :: U_res,U,newU
475
type(scalar_field), intent(inout) :: D_res,D,newD
476
real, dimension(:,:), intent(in) :: local_solver, local_solver_rhs
477
integer, intent(in) :: ele
478
!
479
integer :: uloc, dloc, dim1, d_start, d_end, mdim
480
integer, dimension(mesh_dim(U)) :: U_start, U_end
481
real, dimension(mesh_dim(U)*ele_loc(U,ele)+ele_loc(D,ele))&
482
:: rhs_loc1,rhs_loc2
483
real, dimension(mesh_dim(U),ele_loc(U,ele)) :: U_val
484
485
!Calculate indices in a vector containing all the U and D dofs in
486
!element ele, First the u1 components, then the u2 components, then the
487
!D components are stored.
488
uloc = ele_loc(U_res,ele)
489
dloc = ele_loc(D_res,ele)
490
d_end = mesh_dim(U_res)*uloc+dloc
491
mdim = mesh_dim(U)
492
do dim1 = 1, mdim
493
u_start(dim1) = uloc*(dim1-1)+1
494
u_end(dim1) = uloc*dim1
495
end do
496
d_start = uloc*mdim + 1
497
d_end = uloc*mdim+dloc
498
499
U_val = ele_val(U,ele)
500
do dim1 = 1, mesh_dim(U)
501
rhs_loc1(u_start(dim1):u_end(dim1)) = u_val(dim1,:)
502
end do
503
rhs_loc1(d_start:d_end) = ele_val(D,ele)
504
rhs_loc1 = matmul(local_solver_rhs,rhs_loc1)
505
U_val = ele_val(newU,ele)
506
do dim1 = 1, mesh_dim(U)
507
rhs_loc2(u_start(dim1):u_end(dim1)) = u_val(dim1,:)
508
end do
509
rhs_loc2(d_start:d_end) = ele_val(newD,ele)
510
rhs_loc2 = matmul(local_solver(1:d_end,1:d_end),rhs_loc2)
511
rhs_loc2 = rhs_loc2-rhs_loc1
512
do dim1 = 1, mesh_dim(U)
513
call set(U_res,dim1,ele_nodes(U,ele),&
514
rhs_loc2(u_start(dim1):u_end(dim1)))
515
end do
516
call set(D_res,ele_nodes(D,ele),&
517
rhs_loc2(d_start:d_end))
518
end subroutine get_linear_residuals_ele
519
520
subroutine local_solve_residuals_ele(U_res,D_res,&
521
local_solver_matrix,ele)
522
type(vector_field), intent(inout) :: U_res
523
type(scalar_field), intent(inout) :: D_res
524
real, dimension(:,:), intent(in) :: local_solver_matrix
525
integer, intent(in) :: ele
526
!
527
real, dimension(mesh_dim(U_res)*ele_loc(U_res,ele)+ele_loc(D_res,ele)+&
528
&ele_n_constraints(U_res,ele)) :: rhs_loc
529
real, dimension(mesh_dim(U_res),ele_loc(U_res,ele)) :: U_val
530
integer :: uloc,dloc, d_start, d_end, dim1, mdim
531
integer, dimension(mesh_dim(U_res)) :: U_start, U_end
532
533
rhs_loc = 0.
534
535
!Calculate indices in a vector containing all the U and D dofs in
536
!element ele, First the u1 components, then the u2 components, then the
537
!D components are stored.
538
uloc = ele_loc(U_res,ele)
539
dloc = ele_loc(D_res,ele)
540
d_end = mesh_dim(U_res)*uloc+dloc
541
mdim = mesh_dim(U_res)
542
do dim1 = 1, mdim
543
u_start(dim1) = uloc*(dim1-1)+1
544
u_end(dim1) = uloc*dim1
545
end do
546
d_start = uloc*mdim + 1
547
d_end = uloc*mdim+dloc
548
549
U_val = ele_val(U_res,ele)
550
do dim1 = 1, mesh_dim(U_res)
551
rhs_loc(u_start(dim1):u_end(dim1)) = u_val(dim1,:)
552
end do
553
rhs_loc(d_start:d_end) = &
554
& ele_val(D_res,ele)
555
556
call solve(local_solver_matrix,Rhs_loc)
557
do dim1 = 1, mesh_dim(U_res)
558
call set(U_res,dim1,ele_nodes(U_res,ele),&
559
rhs_loc(u_start(dim1):u_end(dim1)))
560
end do
561
call set(D_res,ele_nodes(D_res,ele),&
562
rhs_loc(d_start:d_end))
563
end subroutine local_solve_residuals_ele
564
565
subroutine assemble_newton_solver_ele(D,f,U,X,down,lambda_rhs,ele, &
566
&g,dt,theta,D0,&
567
&lambda_mat,continuity_mat,&
568
&Local_solver_matrix,Local_solver_rhs)
569
implicit none
570
type(scalar_field), intent(in) :: D,f
571
type(scalar_field), intent(in) :: lambda_rhs
572
type(vector_field), intent(in) :: U,X,down
573
integer, intent(in) :: ele
574
real, intent(in) :: g,dt,theta,D0
575
type(csr_matrix), intent(inout) :: lambda_mat
576
type(block_csr_matrix), intent(inout) :: continuity_mat
577
real, dimension(:,:), intent(inout) ::&
578
& local_solver_matrix
579
real, dimension(:,:), intent(inout) ::&
580
& local_solver_rhs
581
!
582
real, allocatable, dimension(:,:),target :: &
583
&l_continuity_mat, l_continuity_mat2
584
real, allocatable, dimension(:,:) :: helmholtz_loc_mat
585
real, allocatable, dimension(:,:,:) :: continuity_face_mat
586
real, allocatable, dimension(:,:) :: scalar_continuity_face_mat
587
integer :: ni, face
588
integer, dimension(:), pointer :: neigh
589
type(element_type) :: U_shape
590
integer :: stat, d_start, d_end, dim1, mdim, uloc,dloc, lloc, iloc
591
integer, dimension(mesh_dim(U)) :: U_start, U_end
592
type(real_matrix), dimension(mesh_dim(U)) :: &
593
& continuity_mat_u_ptr
594
logical :: have_constraint
595
integer :: constraint_choice, n_constraints, constraints_start
596
597
!Get some sizes
598
lloc = ele_loc(lambda_rhs,ele)
599
mdim = mesh_dim(U)
600
uloc = ele_loc(U,ele)
601
dloc = ele_loc(d,ele)
602
U_shape = ele_shape(U,ele)
603
604
have_constraint = &
605
&have_option(trim(U%mesh%option_path)//"/from_mesh/constraint_type")
606
n_constraints = 0
607
if(have_constraint) then
608
n_constraints = ele_n_constraints(U,ele)
609
end if
610
611
!Calculate indices in a vector containing all the U and D dofs in
612
!element ele, First the u1 components, then the u2 components, then the
613
!D components are stored.
614
do dim1 = 1, mdim
615
u_start(dim1) = uloc*(dim1-1)+1
616
u_end(dim1) = uloc*dim1
617
end do
618
d_start = uloc*mdim + 1
619
d_end = uloc*mdim+dloc
620
621
!Get pointers to different parts of l_continuity_mat
622
allocate(l_continuity_mat(2*uloc+dloc+n_constraints,lloc))
623
do dim1= 1,mdim
624
continuity_mat_u_ptr(dim1)%ptr => &
625
& l_continuity_mat(u_start(dim1):u_end(dim1),:)
626
end do
627
628
! ( M C -L)(u) (v)
629
! ( -C^T N 0 )(h) = (j)
630
! ( L^T 0 0 )(l) (0)
631
!
632
! (u) (M C)^{-1}(v) (M C)^{-1}(L)
633
! (h) = (-C^T N) (j) + (-C^T N) (0)(l)
634
! so
635
! (M C)^{-1}(L) (M C)^{-1}(v)
636
! (L^T 0)(-C^T N) (0)=-(L^T 0)(-C^T N) (j)
637
638
!Get the local_solver matrix that obtains U and D from Lambda on the
639
!boundaries
640
call get_local_solver(local_solver_matrix,U,X,down,D,f,ele,&
641
& g,dt,theta,D0,pullback=.false.,weight=1.0,&
642
& have_constraint=have_constraint,&
643
& projection=.false.,poisson=.false.,&
644
& local_solver_rhs=local_solver_rhs)
645
646
!!!Construct the continuity matrix that multiplies lambda in
647
!!! the U equation
648
!allocate l_continuity_mat
649
l_continuity_mat = 0.
650
!get list of neighbours
651
neigh => ele_neigh(D,ele)
652
!calculate l_continuity_mat
653
do ni = 1, size(neigh)
654
face=ele_face(U, ele, neigh(ni))
655
allocate(continuity_face_mat(mdim,face_loc(U,face)&
656
&,face_loc(lambda_rhs,face)))
657
continuity_face_mat = 0.
658
call get_continuity_face_mat(continuity_face_mat,face,&
659
U,lambda_rhs)
660
do dim1 = 1, mdim
661
continuity_mat_u_ptr(dim1)%ptr(face_local_nodes(U,face),&
662
face_local_nodes(lambda_rhs,face))=&
663
continuity_mat_u_ptr(dim1)%ptr(face_local_nodes(U,face),&
664
face_local_nodes(lambda_rhs,face))+&
665
continuity_face_mat(dim1,:,:)
666
end do
667
do dim1 = 1, mdim
668
call addto(continuity_mat,dim1,1,face_global_nodes(U,face)&
669
&,face_global_nodes(lambda_rhs,face),&
670
&continuity_face_mat(dim1,:,:))
671
end do
672
673
deallocate(continuity_face_mat)
674
end do
675
676
!compute l_continuity_mat2 = inverse(local_solver)*l_continuity_mat
677
allocate(l_continuity_mat2(uloc*2+dloc+n_constraints,lloc))
678
l_continuity_mat2 = l_continuity_mat
679
call solve(local_solver_matrix,l_continuity_mat2)
680
681
!compute helmholtz_loc_mat
682
allocate(helmholtz_loc_mat(lloc,lloc))
683
helmholtz_loc_mat = matmul(transpose(l_continuity_mat),l_continuity_mat2)
684
685
!insert helmholtz_loc_mat into global lambda matrix
686
call addto(lambda_mat,ele_nodes(lambda_rhs,ele),&
687
ele_nodes(lambda_rhs,ele),helmholtz_loc_mat)
688
end subroutine assemble_newton_solver_ele
278
689
279
690
subroutine assemble_hybridized_helmholtz_ele(D,f,U,X,down,ele, &
280
691
g,dt,theta,D0,pullback,weight,lambda_mat,lambda_rhs,U_rhs,D_rhs,&
814
1225
815
1226
end subroutine get_local_solver
816
1227
817
subroutine get_up_gi(X,ele,up_gi,orientation)
818
!subroutine to replace up_gi with a normal to the surface
819
!with the same orientation
820
implicit none
821
type(vector_field), intent(in) :: X
822
integer, intent(in) :: ele
823
real, dimension(X%dim,ele_ngi(X,ele)), intent(inout) :: up_gi
824
integer, intent(out), optional :: orientation
825
!
826
real, dimension(mesh_dim(X), X%dim, ele_ngi(X,ele)) :: J
827
integer :: gi
828
real, dimension(X%dim,ele_ngi(X,ele)) :: normal_gi
829
real, dimension(ele_ngi(X,ele)) :: orientation_gi
830
integer :: l_orientation
831
real :: norm
832
833
call compute_jacobian(ele_val(X,ele), ele_shape(X,ele), J=J)
834
835
select case(mesh_dim(X))
836
case (2)
837
do gi = 1, ele_ngi(X,ele)
838
normal_gi(:,gi) = cross_product(J(1,:,gi),J(2,:,gi))
839
norm = sqrt(sum(normal_gi(:,gi)**2))
840
normal_gi(:,gi) = normal_gi(:,gi)/norm
841
end do
842
do gi = 1, ele_ngi(X,ele)
843
orientation_gi(gi) = dot_product(normal_gi(:,gi),up_gi(:,gi))
844
end do
845
do gi = 1, ele_ngi(X,ele)
846
if(sign(1.0,orientation_gi(gi)).ne.sign(1.0,orientation_gi(1))) then
847
ewrite(0,*) 'gi=',gi
848
ewrite(0,*) 'normal=',normal_gi(:,gi)
849
ewrite(0,*) 'up=',up_gi(:,gi)
850
ewrite(0,*) 'orientation=',orientation_gi(gi)
851
ewrite(0,*) 'orientation(1)=',orientation_gi(1)
852
FLAbort('Nasty geometry problem')
853
end if
854
end do
855
856
857
if(orientation_gi(1)>0.0) then
858
l_orientation = 1
859
else
860
l_orientation = -1
861
end if
862
if(present(orientation)) then
863
orientation =l_orientation
864
end if
865
do gi = 1, ele_ngi(X,ele)
866
up_gi(:,gi) = normal_gi(:,gi)*l_orientation
867
end do
868
case default
869
FLAbort('not implemented')
870
end select
871
end subroutine get_up_gi
872
873
1228
function get_orientation(X_val, up) result (orientation)
874
1229
!function compares the orientation of the element with the
875
1230
!up direction
1417
1772
type(vector_field), pointer :: U_local,down,X, U_cart
1418
1773
type(vector_field) :: Coriolis_term, Balance_eqn, tmpV_field
1419
1774
integer :: ele,dim1,i1, stat
1420
real :: g
1775
real :: g, D0
1421
1776
logical :: elliptic_method, pullback
1422
1777
real :: u_max, b_val, h_mean, area
1423
1778
real, dimension(:), allocatable :: weights
1514
1869
end do
1515
1870
h_mean = h_mean/area
1516
1871
D%val = D%val - h_mean
1872
!Add back on the correct mean depth
1873
call get_option("/material_phase::Fluid/scalar_field::LayerThickness/&
1874
&prognostic/mean_layer_thickness",D0)
1875
D%val = D%val + D0
1517
1876
1518
1877
!debugging tests
1519
1878
call zero(Coriolis_term)
1759
2118
real, intent(in) :: weight
1760
2119
!
1761
2120
real, dimension(ele_loc(psi,ele)) :: psi_loc
1762
real, dimension(mesh_dim(psi),ele_ngi(psi,ele)) :: dpsi_gi
2121
real, dimension(mesh_dim(psi),ele_ngi(psi,ele)) :: dpsi_gi, Uloc_gi
1763
2122
real, dimension(ele_ngi(psi,ele)) :: div_gi
1764
2123
real, dimension(mesh_dim(U_local),ele_loc(U_local,ele)) :: U_loc
1765
2124
real, dimension(ele_loc(U_local,ele),ele_loc(U_local,ele)) :: &
1801
2160
call solve(l_mass_mat,U_loc(dim1,:))
1802
2161
end do
1803
2162
2163
!DEbugging check, did we get the same fields?
2164
do dim1= 1, U_local%dim
2165
Uloc_gi(dim1,:) = matmul(transpose(U_shape%n),U_loc(dim1,:))
2166
end do
2167
assert(maxval(abs(Uloc_gi-dpsi_gi))/max(1.0,maxval(abs(Uloc_gi)))<1.0e-8)
2168
1804
2169
!verify divergence-free-ness
1805
2170
!This is just for debugging
1806
2171
!Annoyingly it requires detJ when the rest of the subroutine doesn't
1894
2259
D_gi = ele_val_at_quad(D,ele)
1895
2260
call compute_jacobian(ele_val(X,ele),ele_shape(X,ele), J=J, &
1896
2261
detwei=detwei)
1897
2262
assert(all(detwei>0.))
1898
2263
mean = mean + sum(detwei*D_gi)
1899
2264
area = area + sum(detwei)
1900
2265
2131
2496
end if
2132
2497
2133
2498
!construct a tangent basis on the manifold
2499
FLExit('This is bugged because vectors arent normalised.')
2134
2500
do gi = 1, ele_ngi(X,ele)
2135
2501
m_basis(1,:,gi) = cross_product(ref_vec_gi(:,gi),up_gi(:,gi))
2136
2502
m_basis(2,:,gi) = cross_product(up_gi(:,gi),m_basis(1,:,gi))
2324
2690
2325
2691
if(present(name)) then
2326
2692
D=>extract_scalar_field(state, trim(name))
2327
call get_option("/material_phase::Fluid/scalar_field::PrescribedLayer&
2328
&DepthFromProjection/prescribed/python",Python_Function)
2693
call get_option(trim(D%option_path)//"/prescribed/value/python",&
2694
& Python_Function)
2329
2695
else
2330
2696
D=>extract_scalar_field(state, "LayerThickness")
2331
2697
call get_option("/material_phase::Fluid/scalar_field&
2339
2705
call set_layerthickness_projection_ele(D,X,Python_Function,ele)
2340
2706
end do
2341
2707
2342
!Subtract off the mean part
2343
h_mean = 0.
2344
area = 0.
2345
do ele = 1, element_count(D)
2346
call assemble_mean_ele(D,X,h_mean,area,ele)
2347
end do
2348
h_mean = h_mean/area
2349
D%val = D%val - h_mean
2350
2351
2708
end subroutine set_layerthickness_projection
2352
2709
2353
2710
subroutine set_layerthickness_projection_ele(D,X,Python_Function,ele)
2479
2836
2480
2837
end subroutine set_velocity_from_lat_long_ele
2481
2838
2482
subroutine set_velocity_from_sphere_pullback(state)
2839
subroutine set_velocity_from_sphere_pullback(state,v_field&
2840
&,Python_Function_in,R0_in)
2483
2841
type(state_type), intent(inout) :: state
2842
type(vector_field), target, intent(inout), optional :: v_field
2843
character(len=*), intent(in), optional :: Python_Function_in
2844
real, intent(in), optional :: R0_in
2484
2845
!
2485
2846
type(vector_field), pointer :: X,U
2486
2847
character(len=PYTHON_FUNC_LEN) :: Python_Function
2488
2849
integer :: ele
2489
2850
2490
2851
X=>extract_vector_field(state, "Coordinate")
2491
U=>extract_vector_field(state, "LocalVelocity")
2492
call get_option("/material_phase::Fluid/vector_field::Velocity&
2493
&/prognostic/initial_condition::WholeMesh/&
2494
&from_sphere_pullback/python",Python_Function)
2495
call get_option("/material_phase::Fluid/vector_field::Velocity&
2496
&/prognostic/initial_condition::WholeMesh/&
2497
&from_sphere_pullback/sphere_radius",R0)
2852
if(present(v_field)) then
2853
U=>v_field
2854
else
2855
U=>extract_vector_field(state, "LocalVelocity")
2856
end if
2857
if(present(Python_Function_in)) then
2858
Python_Function = trim(Python_Function_in)
2859
else
2860
call get_option("/material_phase::Fluid/vector_field::Velocity&
2861
&/prognostic/initial_condition::WholeMesh/&
2862
&from_sphere_pullback/python",Python_Function)
2863
end if
2864
if(present(R0_in)) then
2865
R0 = R0_in
2866
else
2867
call get_option("/material_phase::Fluid/vector_field::Velocity&
2868
&/prognostic/initial_condition::WholeMesh/&
2869
&from_sphere_pullback/sphere_radius",R0)
2870
end if
2498
2871
2499
2872
do ele = 1, element_count(U)
2500
2873
call set_velocity_from_sphere_pullback_ele&
2558
2931
m_normal(1,:) = xm/(xm**2+ym**2+zm**2)**0.5
2559
2932
m_normal(2,:) = ym/(xm**2+ym**2+zm**2)**0.5
2560
2933
m_normal(3,:) = zm/(xm**2+ym**2+zm**2)**0.5
2561
assert(maxval(abs(sum(Us_quad*m_normal,1)))<1.0e-7)
2562
2934
assert(maxval(abs(sum(m_normal**2,1)-1.0))<1.0e-8)
2935
if(maxval(abs(sum(Us_quad*m_normal,1)))>1.0e-7) then
2936
do gi = 1, ele_ngi(X,ele)
2937
ewrite(2,*) 'Us',Us_quad(:,gi)
2938
ewrite(2,*) 'normal', m_normal(:,gi)
2939
ewrite(2,*) 'dot prod', sum(Us_quad(:,gi)*m_normal(:,gi))
2940
end do
2941
FLExit('Velocity field not tangential to sphere')
2942
end if
2563
2943
e_normal(1,:) = xe/(xe**2+ye**2+ze**2)**0.5
2564
2944
e_normal(2,:) = ye/(xe**2+ye**2+ze**2)**0.5
2565
2945
e_normal(3,:) = ze/(xe**2+ye**2+ze**2)**0.5
2658
3038
do dim1 = 1, mesh_dim(U)
2659
3039
U_local_loc(dim1,:) = l_u_rhs((dim1-1)*uloc+1:dim1*uloc)
2660
3040
end do
2661
!debugging stuff
2662
U_local_quad = matmul(U_local_loc,u_shape%n)
2663
do gi = 1, ele_ngi(U,ele)
2664
U_quad_2(:,gi) = matmul(transpose(J(:,:,gi)),U_local_quad(:,gi))/detJ(gi)
2665
end do
2666
assert(maxval(abs(U_quad-U_quad_2))<1.0e-7)
2667
3041
2668
3042
do dim1 = 1, mesh_dim(U)
2669
3043
call set(U,ele_nodes(U,ele),U_local_loc)
Loggerhead is a web-based interface for Breezy
Version: 2.0.1