« back to all changes in this revision
Viewing changes to assemble/Momentum_Equation.F90
- browse files at revision 3999
- compare with another revision
- download diff
- download tarball
- view history from revision 3999
- Committer: traverscug
- Date: 2012-11-12 18:41:00 UTC
- Revision ID: traverscug@gmail.com-20121112184100-os0j9q944e2ezarr
Uncommit last revisions
- files modified:
- adjoint_reduced_modelling/Adjoint_main.F90
added
removed
assemble/Momentum_Equation.F90
122
122
logical :: diagonal_big_m
123
123
logical :: pressure_debugging_vtus
124
124
!! True if the momentum equation should be solved with the reduced model.
125
logical :: reduced_model,petrov
125
logical :: reduced_model
126
126
127
127
! Add viscous terms to inverse_masslump for low Re which is only used for pressure correction
128
128
logical :: low_re_p_correction_fix
263
263
! Do we have fluid-particle drag between phases?
264
264
logical :: have_fp_drag
265
265
266
267
!petro
268
type(scalar_field), pointer :: POD_u_scalar !petro
269
type(pod_rhs_type)::pod_rhs_old
270
type(vector_field), pointer :: POD_u
271
real, dimension(:), allocatable :: THETA,smean_gi,KB_pod,b_pod,PSI_GI,PSI_OLD_GI,psi_old,GRADX,GRADT,DIFFGI_pod,PSI
272
real, dimension(:,:), allocatable ::leftsvd,leftsvd_gi,leftsvd_x_gi, AMAT_pod, KMAT_pod
273
274
!real, dimension(:,:), allocatable :: snapmean_velocity
275
type(element_type), pointer :: x_shape,xf_shape
276
real, dimension(:,:),allocatable ::X_val !for getting the result of detwei
277
real, dimension(:), allocatable :: detwei(:),res(:)
278
real ,dimension(:,:),allocatable :: JJ !(X%dim, mesh_dim(X)),
279
real :: det ! for getting the result of detwei
280
real, dimension(:,:,:), allocatable :: dshape
281
282
! real, dimension(POD_u%dim,ele_loc(POD_u,ele)) :: X_val !for getting the result of detwei
283
! new added for petro-galerkin
284
integer :: TOTELE,NLOC,NGI,nsvd !,stat
285
integer :: POD_num
286
type(mesh_type), pointer :: pod_umesh
287
REAL PSI_THETA
288
INTEGER ELE,globi,globj,isvd,jsvd,jloc,GI
289
REAL A_STAR_COEF,AX_STAR,P_STAR_POD
290
REAL DT1
291
real noloc,nogas
292
293
294
266
295
267
!!for reduced model
296
268
type(vector_field), pointer :: snapmean_velocity
297
269
type(scalar_field), pointer :: snapmean_pressure
329
301
330
302
logical :: on_sphere, have_absorption, have_vertical_stabilization
331
303
type(vector_field), pointer :: dummy_absorption
332
304
333
305
ewrite(1,*) 'Entering solve_momentum'
334
306
call get_option('/timestepping/nonlinear_iterations', nonlinear_iterations, default=1)
335
307
340
312
endif
341
313
342
314
! reduced_model= have_option("/reduced_model/execute_reduced_model")
343
petrov=.true.
344
DT1=0.1
315
345
316
! if(timestep==1)timestep_check=.true.
346
317
if(timestep==1)then
347
318
timestep_check=.true.
1326
1297
enddo
1327
1298
1328
1299
!100 continue
1329
!test
1330
if(petrov) then
1331
print *,"before POD_u"
1332
POD_u=>extract_vector_field(POD_state(1,1,1), "Velocity")
1333
print *,"before POD_u_scalar"
1334
!POD_u_scalar=extract_scalar_field(POD_u,1) !!need modify
1335
POD_num=size(POD_state,1)
1336
print *,"before allocate(X_val"
1337
allocate(X_val(POD_u%dim,ele_loc(POD_u,ele)))
1338
print *,"before allocate JJ"
1339
allocate(JJ(POD_u%dim,mesh_dim(POD_u)))
1340
!allocate(pod_coef(POD_u%dim*POD_num+POD_num))
1341
pod_umesh => extract_mesh(pod_state(1,1,1), "Mesh", stat)
1342
TOTELE=pod_umesh%elements
1343
NLOC=pod_umesh%shape%loc
1344
NGI=pod_umesh%shape%ngi
1345
print *,"before allocate rhs_old"
1346
call allocate(pod_rhs_old, POD_u, POD_pressure)
1347
1348
xf_shape=>ele_shape(POD_u,1)
1349
x_shape=>ele_shape(POD_u,1)
1350
print *, "timestep,, 2+++"
1351
do j=1, POD_num
1352
!POD_u_scalar=>extract_scalar_field(POD_state(j), "PODPressure")
1353
!POD_u_scalar=extract_scalar_field(POD_u,j)
1354
X_val=ele_val(POD_u, j)
1355
xf_shape=>ele_shape(Pod_u,j) !!
1356
x_shape=>ele_shape(POD_u,j)!!
1357
end do
1358
!!new added for petro-galerkin
1359
1360
1361
! allocate(pod_sol_velocity(u_nodes,POD_velocity%dim))
1362
allocate(THETA(NGI))
1363
allocate(leftsvd(TOTELE*NLOC+NLOC,nsvd))
1364
allocate(leftsvd_gi(NGI,nsvd))
1365
allocate(leftsvd_x_gi(NGI,nsvd))
1366
allocate(smean_gi(NGI))
1367
allocate(AMAT_pod(nsvd,nsvd))
1368
allocate(KMAT_pod(nsvd,nsvd))
1369
allocate(KB_pod(nsvd))
1370
allocate(b_pod(nsvd))
1371
allocate(PSI(TOTELE*NLOC+NLOC))
1372
allocate(PSI_GI(NGI))
1373
allocate(PSI_OLD_GI(NGI))
1374
allocate(psi_old(TOTELE*NLOC+NLOC))
1375
allocate(GRADX(NGI))
1376
allocate(GRADT(NGI))
1377
allocate(DIFFGI_pod(NGI))
1378
allocate(detwei(NGI))
1379
allocate(res(NGI))
1380
! noloc=x_shape%loc
1381
! nogas=x_shape%ngi
1382
! allocate(dshape(noloc,nogas,POD_u%dim))
1383
allocate(dshape(x_shape%loc,x_shape%ngi,POD_u%dim)) !number of nodes, NO. of gussian
1384
1385
1386
! new added for petro-galerkin Residual, see Non-Linear Petrov-Galerkin Methods for
1387
! Hyperbolic Equations and Discontinuous Finite Element Methods .. eqn 102
1388
! get the result of Residual value"
1389
1390
pod_rhs_old%val=pod_rhs%val
1391
call solve(pod_matrix%val, pod_rhs%val)
1392
! pod_coef=pod_rhs%val
1393
pod_rhs%val=pod_rhs%val-pod_rhs_old%val
1394
call SMLINN(pod_matrix%val,RES,pod_rhs%val,GI,GI) ! get the RES
1395
1396
DO ELE=1,TOTELE
1397
leftsvd_gi = 0.0
1398
leftsvd_x_gi = 0.0
1399
smean_gi =0.0
1400
1401
call transform_to_physical(POD_u, ele, ele_shape(POD_u,1), detwei=detwei, dshape=dshape)
1402
DO GI=1,NGI ! guassian points' number
1403
DO JLOC=1,NLOC ! NLOC points in one element
1404
GLOBJ=(ELE-1)*NLOC+JLOC
1405
1406
do isvd =1,nsvd
1407
leftsvd(GLOBJ,isvd)=POD_u%val(GLOBJ,isvd)
1408
leftsvd_gi(GI,isvd)=leftsvd_gi(GI,isvd)+xf_shape%dn(isvd,GI,1)*leftsvd(GLOBJ,isvd)
1409
!N shape function, leftsvd: podbasis, N(JLOC,GI): ele_shape(nu, ele)
1410
leftsvd_x_gi(GI,isvd)=leftsvd_gi(GI,isvd)+dshape(JLOC,GI,1)*leftsvd(GLOBJ,isvd)
1411
enddo
1412
!smean_gi(gi)=smean_gi(gi)+ 1SMLINN(A,X,B,NMX,N)*smean(GLOBJ) !need modification
1413
!smean_gi(gi)=smean_gi(gi)+ 1*smean(GLOBJ)
1414
END DO
1415
END DO
1416
END DO
1417
1418
PSI_GI=0.0 ! quadrature point
1419
PSI_OLD_GI=0.0
1420
1421
DO GI=1,NGI
1422
DO JLOC=1,NLOC
1423
GLOBJ=(ELE-1)*NLOC+JLOC
1424
PSI_GI(GI)=PSI_GI(GI)+xf_shape%dn(isvd,GI,1)*PSI(GLOBJ)
1425
PSI_OLD_GI(GI)=PSI_OLD_GI(GI)+xf_shape%dn(isvd,GI,1)*PSI_OLD(GLOBJ)
1426
END DO
1427
END DO
1428
1429
1430
! calculate the GRADX and GRADT quadrature point
1431
GRADX=0.0
1432
GRADT=0.0
1433
DO GI=1,NGI ! ?? NGI and NLOC need input
1434
DO JLOC=1,NLOC
1435
GLOBJ=(ELE-1)*NLOC+JLOC
1436
PSI_THETA=THETA(GI)*PSI(GLOBJ)+(1.-THETA(GI))*PSI_OLD(GLOBJ)
1437
1438
GRADX(GI)=GRADX(GI)+dshape(JLOC,GI,1)*PSI_THETA
1439
GRADT(GI)=GRADT(GI)+xf_shape%dn(isvd,GI,1)*(PSI(GLOBJ)-PSI_OLD(GLOBJ))/DT1
1440
END DO
1441
END DO
1442
1443
! get the result of AX_STAR cos P_STAR_POD need it
1444
DO GI=1,NGI
1445
1446
A_STAR_COEF=(GRADX(GI)+GRADT(GI)*1.) &
1447
/MAX(1.E-6,GRADX(GI)**2+1.*GRADT(GI)**2)
1448
! AT_STAR=A_STAR_COEF*GRADT(GI)
1449
AX_STAR=A_STAR_COEF*GRADX(GI)
1450
enddo
1451
do gi =1,ngi
1452
do isvd = 1,nsvd
1453
if(isvd.eq.1) then
1454
P_STAR_POD = 0.25/abs(AX_STAR*leftsvd_x_gi(GI,isvd))
1455
else
1456
P_STAR_POD = 0.25/max(P_STAR_POD,abs(AX_STAR*leftsvd_x_gi(GI,isvd)))
1457
endif
1458
enddo
1459
DIFFGI_pod(GI) = 1.0e-4*P_STAR_POD*RES(GI)**2 &
1460
/ MAX(1.E-6,GRADX(GI)**2+GRADT(GI)**2)!!v_pod (eqn.100)
1461
enddo
1462
1463
do isvd =1,nsvd
1464
do jsvd=1,nsvd
1465
do gi =1,ngi
1466
KMAT_pod(isvd,jsvd) = KMAT_pod(isvd,jsvd) &
1467
+ detwei(GI)*leftsvd_x_gi(GI,isvd)*DIFFGI_pod(GI)*leftsvd_x_gi(GI,jsvd)
1468
! (eqn.99 left side,leftsvd_x_gi:delta N)
1469
enddo
1470
enddo
1471
enddo
1472
1473
1474
!!call solver
1475
1476
1477
1478
!pod_rhs%val=0.001*pod_rhs%val
1479
! add K to (MT AM + K) fi *pod= MT b.
1480
1481
do isvd =1,nsvd
1482
do jsvd=1,nsvd
1483
pod_matrix%val(isvd,jsvd) = pod_matrix%val(isvd,jsvd)+KMAT_pod(isvd,jsvd)
1484
enddo
1485
!b_pod(isvd) = b_pod(isvd)-KB_pod(isvd)
1486
enddo
1487
1488
call solve(pod_matrix%val,pod_rhs_old%val)
1489
pod_rhs%val=pod_rhs_old%val
1490
else
1491
call solve(pod_matrix%val,pod_rhs%val)
1492
! print *, 'size(pod_rhs%val)size(pod_rhs%val)size(pod_rhs%val)', size(pod_rhs%val)
1493
endif ! end petrov
1494
! call solve(pod_matrix%val, pod_rhs%val)
1300
1301
call solve(pod_matrix%val, pod_rhs%val)
1495
1302
pod_coef_dt=0.0
1496
1303
pod_coef_dt=pod_rhs%val
1497
1304
pod_rhs%val=0.0
2573
2380
ewrite(1,*) 'Finished checking momentum discretisation options'
2574
2381
2575
2382
end subroutine momentum_equation_check_options
2576
SUBROUTINE SMLINN(A,X,B,NMX,N)
2577
IMPLICIT NONE
2578
INTEGER NMX,N
2579
REAL A(NMX,NMX),X(NMX),B(NMX)
2580
REAL R
2581
INTEGER K,I,J
2582
! Form X = A^{-1} B
2583
! Useful subroutine for inverse
2584
! This sub overwrites the matrix A.
2585
DO K=1,N-1
2586
DO I=K+1,N
2587
A(I,K)=A(I,K)/A(K,K)
2588
END DO
2589
DO J=K+1,N
2590
DO I=K+1,N
2591
A(I,J)=A(I,J) - A(I,K)*A(K,J)
2592
END DO
2593
END DO
2594
END DO
2595
!
2596
! Solve L_1 x=b
2597
DO I=1,N
2598
R=0.
2599
DO J=1,I-1
2600
R=R+A(I,J)*X(J)
2601
END DO
2602
X(I)=B(I)-R
2603
END DO
2604
!
2605
! Solve U x=y
2606
DO I=N,1,-1
2607
R=0.
2608
DO J=I+1,N
2609
R=R+A(I,J)*X(J)
2610
END DO
2611
X(I)=(X(I)-R)/A(I,I)
2612
END DO
2613
RETURN
2614
END SUBROUTINE SMLINN
2383
2615
2384
end module momentum_equation
Loggerhead is a web-based interface for Breezy
Version: 2.0.1