← Back to branch summary

~reducedmodelling/fluidity/ROM_DA_DEIM

« back to all changes in this revision

Viewing changes to assemble/Momentum_Equation.F90

  • Committer: fangf at ac
  • Date: 2013-01-10 13:28:22 UTC
  • Revision ID: fangf@imperial.ac.uk-20130110132822-s7l0zb9n0mkhhnea
tidy up

Show diffs side-by-side

added added

removed removed

Lines of Context:
assemble/Momentum_Equation.F90
687
687
            
688
688
            ! If CV pressure then add in any dirichlet pressure BC integrals to the mom_rhs.
689
689
            if (cv_pressure) then
 
690
!            if (cv_pressure.and.(.not.reduced_model)) then
690
691
               call add_pressure_dirichlet_bcs_cv(mom_rhs(istate), u, p, state(istate))
691
692
            end if
692
693
            
1168
1169
 
1169
1170
            call get_option('/timestepping/nonlinear_iterations', nonlinear_iterations, default=1)
1170
1171
            !!print*,'its = ', its, ' nonlinear_iterations = ', nonlinear_iterations
1171
 
 
 
1172
            
1172
1173
            !enter reduced model process after the last nonliear_iteration
1173
 
           ! if(its==nonlinear_iterations)then
1174
 
               call profiler_tic("reduced_model_loop")
1175
 
 
1176
 
               reduced_model_loop: do istate = 1, size(state)
1177
 
 
 
1174
            ! if(its==nonlinear_iterations)then
 
1175
            call profiler_tic("reduced_model_loop")
 
1176
            
 
1177
            reduced_model_loop: do istate = 1, size(state)
1178
1178
               ! Get the velocity
1179
1179
               u => extract_vector_field(state(istate), "Velocity", stat)
1180
1180
               lump_mass_form = have_option(trim(u%option_path)//&
1190
1190
               if(.not.have_option(trim(u%option_path)//"/prognostic")) cycle
1191
1191
 
1192
1192
               if(have_option(trim(u%option_path)//"/prognostic/spatial_discretisation"//&
1193
 
                                       &"/continuous_galerkin").or.&
1194
 
                  have_option(trim(u%option_path)//"/prognostic/spatial_discretisation"//&
1195
 
                                       &"/discontinuous_galerkin")) then
1196
 
 
1197
 
!!speedup version of reduced model
1198
 
 
 
1193
                    &"/continuous_galerkin").or.&
 
1194
                    have_option(trim(u%option_path)//"/prognostic/spatial_discretisation"//&
 
1195
                    &"/discontinuous_galerkin")) then
 
1196
                  
 
1197
                  !!speedup version of reduced model
 
1198
                  
1199
1199
                  ! Allocate the change in pressure field
1200
1200
                  call allocate(delta_p, p%mesh, "DeltaP")
1201
1201
                  delta_p%option_path = trim(p%option_path)
1202
1202
                  call zero(delta_p)
1203
 
 
 
1203
                  
1204
1204
                  call allocate(delta_u, u%dim, u%mesh, "DeltaU")
1205
1205
                  delta_u%option_path = trim(u%option_path)
1206
1206
                  call zero(delta_u)
1207
 
 
1208
 
 
 
1207
                  
 
1208
                  
1209
1209
                  call get_option("/timestepping/timestep", dt)
1210
 
               
 
1210
                  
1211
1211
                  POD_velocity=>extract_vector_field(POD_state(1,1,istate), "Velocity")
1212
1212
                  POD_pressure=>extract_scalar_field(POD_state(1,2,istate), "Pressure")    
1213
1213
                  call allocate(pod_matrix, POD_velocity, POD_pressure)
1214
1214
                  call allocate(pod_matrix_mass, POD_velocity, POD_pressure)
1215
1215
                  call allocate(pod_matrix_adv, POD_velocity, POD_pressure)             
1216
1216
                  call allocate(pod_rhs, POD_velocity, POD_pressure)
1217
 
                   
 
1217
                  
1218
1218
                  allocate(pod_matrix_snapmean((u%dim+1)*size(POD_state,1),(u%dim+1)*size(POD_state,1)))
1219
1219
                  allocate(pod_matrix_adv_mean((u%dim+1)*size(POD_state,1),(u%dim+1)*size(POD_state,1)))
1220
1220
                  allocate(pod_rhs_snapmean((u%dim+1)*size(POD_state,1)))
1228
1228
 
1229
1229
                
1230
1230
                  if (timestep==1)then
1231
 
                        !!the ct_rhs for reduced model
1232
 
                         ct_rhs(istate)%val=ct_rhs(istate)%val/dt/theta_divergence
1233
 
 
1234
 
                         ! print*,'theta_pg',theta_pg
1235
 
                         ! print*,'theta_div',theta_divergence
1236
 
                         ! print*,'ct_rhs',ct_rhs(istate)%val(1:10)
1237
 
                          
1238
 
 
 
1231
                     !!the ct_rhs for reduced model
 
1232
                     ct_rhs(istate)%val=ct_rhs(istate)%val/dt/theta_divergence
 
1233
                     
 
1234
                     ! print*,'theta_pg',theta_pg
 
1235
                     ! print*,'theta_div',theta_divergence
 
1236
                     ! print*,'ct_rhs',ct_rhs(istate)%val(1:10)
 
1237
                     
 
1238
                     
1239
1239
                     print*,has_boundary_condition(u,'free_surface')
1240
1240
                     !project to reduced matrix
1241
1241
                     if(has_boundary_condition(u,'free_surface'))then
1242
 
                         call project_reduced(big_m(istate), mom_rhs(istate), ct_m(istate)%ptr, ct_rhs(istate), pod_matrix, &
 
1242
                        call project_reduced(big_m(istate), mom_rhs(istate), ct_m(istate)%ptr, ct_rhs(istate), pod_matrix, &
1243
1243
                             pod_rhs, dt, POD_state(:,:,istate), fs_m=fs_m)
1244
1244
                     else
1245
 
                         call project_reduced(big_m(istate), mom_rhs(istate), ct_m(istate)%ptr, ct_rhs(istate), pod_matrix, &
 
1245
                        call project_reduced(big_m(istate), mom_rhs(istate), ct_m(istate)%ptr, ct_rhs(istate), pod_matrix, &
1246
1246
                             pod_rhs, dt, POD_state(:,:,istate))
1247
1247
                     endif
1248
 
 
 
1248
                     
1249
1249
                     !save the pod_matrix and pod_rhs for snapmean state
1250
1250
                     if(snapmean)then
1251
1251
                        ewrite(1,*)'snapmean'
1253
1253
                        write(20,*)((pod_matrix%val(i,j),j=1,(u%dim+1)*size(POD_state,1)),i=1,(u%dim+1)*size(POD_state,1))
1254
1254
                        write(20,*)(pod_rhs%val(i),i=1,(u%dim+1)*size(POD_state,1))
1255
1255
                        close(20)
1256
 
                                  
1257
 
                         call zero(big_m_tmp(istate))
 
1256
                        
 
1257
                        call zero(big_m_tmp(istate))
1258
1258
                        pod_rhs%val=0.0
1259
 
                      !  pod_matrix%val=0.0
 
1259
                        !  pod_matrix%val=0.0
1260
1260
                        pod_matrix_mass%val=0.0
1261
1261
                        pod_matrix_adv%val=0.0
1262
1262
                        if(lump_mass_form) then
1267
1267
                                 !! we comment invert(inverse_masslump) in Momentum_CG
1268
1268
                                 !! here, big_m_tmp is tmp matrix
1269
1269
                                 call addto_diag(big_m_tmp(istate), dim, dim,i,inverse_masslump(istate)%val(i,dim) )
1270
 
                             !  print *,'inverse_masslump(istate)%val(i,dim)', inverse_masslump(istate)%val(i,dim)
1271
 
                                
 
1270
                                 !  print *,'inverse_masslump(istate)%val(i,dim)', inverse_masslump(istate)%val(i,dim)
 
1271
                                 
1272
1272
                              enddo
1273
1273
                           enddo
1274
 
                          call project_reduced(big_m_tmp(istate), mom_rhs(istate), ct_m(istate)%ptr, ct_rhs(istate), pod_matrix_mass, &
1275
 
                             pod_rhs, dt, POD_state(:,:,istate))
 
1274
                           call project_reduced(big_m_tmp(istate), mom_rhs(istate), ct_m(istate)%ptr, ct_rhs(istate), pod_matrix_mass, &
 
1275
                                pod_rhs, dt, POD_state(:,:,istate))
1276
1276
                           !  print *,'pod_matrix_mass', pod_matrix_mass%val
1277
 
                             
 
1277
                           
1278
1278
                        else
1279
 
                          ! print*,'mass'
 
1279
                           ! print*,'mass'
1280
1280
                           !call project_reduced(mass(istate), mom_rhs(istate), ct_m(istate)%ptr, ct_rhs(istate), pod_matrix, &
1281
 
                             !   pod_ct_m, pod_rhs, 1.0, POD_state(:,:,istate))
 
1281
                           !   pod_ct_m, pod_rhs, 1.0, POD_state(:,:,istate))
1282
1282
                        endif
1283
 
                       
1284
 
                         !save the mass_matrix 
 
1283
                        
 
1284
                        !save the mass_matrix 
1285
1285
                        open(unit=20,file='mass_matrix')
1286
1286
                        write(20,*)((pod_matrix_mass%val(i,j),j=1,(u%dim+1)*size(POD_state,1)),i=1,(u%dim+1)*size(POD_state,1))
1287
1287
                        close(20)
1289
1289
                        !save the advection_matrix 
1290
1290
                        pod_matrix_adv%val(:,:)=pod_matrix%val(:,:)-pod_matrix_mass%val(:,:)
1291
1291
                        open(unit=21,file='advection_matrix_snapmean')
1292
 
                         write(21,*)((pod_matrix_adv%val(i,j),j=1,(u%dim+1)*size(POD_state,1)),i=1,(u%dim+1)*size(POD_state,1))
1293
 
                         close(21)
1294
 
                       !  open(unit=21,file='big_m')
1295
 
                       !  write(21,*)((pod_matrix%val(i,j),j=1,(u%dim+1)*size(POD_state,1)),i=1,(u%dim+1)*size(POD_state,1))
1296
 
                      !   close(21)
1297
 
        
1298
 
                      !  open(unit=20,file='ct_m_matrix')
1299
 
                       ! write(20,*) ((pod_ct_m(i,j),j=1,u%dim*size(POD_state,1)),i=1,size(POD_state,1))
1300
 
                      !  close(20)
1301
 
                     endif      
 
1292
                        write(21,*)((pod_matrix_adv%val(i,j),j=1,(u%dim+1)*size(POD_state,1)),i=1,(u%dim+1)*size(POD_state,1))
 
1293
                        close(21)
 
1294
                        !  open(unit=21,file='big_m')
 
1295
                        !  write(21,*)((pod_matrix%val(i,j),j=1,(u%dim+1)*size(POD_state,1)),i=1,(u%dim+1)*size(POD_state,1))
 
1296
                        !   close(21)
 
1297
                        
 
1298
                        !  open(unit=20,file='ct_m_matrix')
 
1299
                        ! write(20,*) ((pod_ct_m(i,j),j=1,u%dim*size(POD_state,1)),i=1,size(POD_state,1))
 
1300
                        !  close(20)
 
1301
                     endif
1302
1302
                     
1303
1303
                     if(eps.ne.0.0)then
1304
1304
                        pod_matrix_mass%val=0.0
1315
1315
                        pod_matrix_adv_perturbed(:,:)=(pod_matrix_adv%val(:,:)-pod_matrix_adv_mean(:,:))/eps
1316
1316
                        !save the advection_matrix 
1317
1317
                        ! open(unit=21,file='advection_matrix_perturbed')
1318
 
                         write(60,*)((pod_matrix_adv_perturbed(i,j),j=1,(u%dim+1)*size(POD_state,1)),i=1,(u%dim+1)*size(POD_state,1))
1319
 
                       !  close(21)
1320
 
                         
 
1318
                        write(60,*)((pod_matrix_adv_perturbed(i,j),j=1,(u%dim+1)*size(POD_state,1)),i=1,(u%dim+1)*size(POD_state,1))
 
1319
                        !  close(21)
 
1320
                        
1321
1321
                        ewrite(1,*)'perturbed'
1322
1322
                        print*,'perturbed'
1323
1323
                        open(unit=20,file='pod_matrix_snapmean')
1324
1324
                        read(20,*)((pod_matrix_snapmean(i,j),j=1,(u%dim+1)*size(POD_state,1)),i=1,(u%dim+1)*size(POD_state,1))
1325
1325
                        read(20,*)(pod_rhs_snapmean(i),i=1,(u%dim+1)*size(POD_state,1))
1326
1326
                        close(20)
1327
 
 
1328
 
                        pod_matrix%val(:,:)=(pod_matrix%val(:,:)-pod_matrix_snapmean(:,:))/eps
 
1327
                        pod_matrix%val(:,:)=(pod_matrix%val(:,:)-pod_matrix_snapmean(:,:))/eps
1329
1328
                        pod_rhs%val(:)=(pod_rhs%val(:)-pod_rhs_snapmean(:))/eps
1330
 
                      
 
1329
                        
1331
1330
                        do i=1,size(POD_state,1)
1332
1331
                           POD_velocity=>extract_vector_field(POD_state(i,1,istate), "Velocity")
1333
1332
                           do j=1,u%dim
1336
1335
                           POD_pressure=>extract_scalar_field(POD_state(i,2,istate), "Pressure")
1337
1336
                           pod_coef(i+size(POD_state,1)*u%dim)=dot_product(POD_pressure,p)
1338
1337
                        enddo
1339
 
                        call Matrix_vector_multiplication(size(pod_coef),pod_rhs_adv_perturbed, &
1340
 
                                                pod_matrix_adv_perturbed,pod_coef)
 
1338
                        call Matrix_vector_multiplication(size(pod_coef),pod_rhs_adv_perturbed, &
 
1339
                             pod_matrix_adv_perturbed,pod_coef)
1341
1340
                        
1342
1341
                        pod_rhs%val = pod_rhs%val + pod_rhs_adv_perturbed ! delete advection term in rhs
1343
 
 
 
1342
   
1344
1343
                        !save pod_matrix and pod_rhs for perturbed state (related to POD_velocity%dim, size(POD_state,1))
1345
1344
                        write(30,*)((pod_matrix%val(i,j),j=1,(u%dim+1)*size(POD_state,1)),i=1,(u%dim+1)*size(POD_state,1))
1346
1345
                        write(50,*)(pod_rhs%val(i),i=1,(u%dim+1)*size(POD_state,1))! delete advection term in rhs
1348
1347
                     elseif(eps.eq.0.0 .and. .not.snapmean)then
1349
1348
                        !need to exclude snapmean
1350
1349
                        ewrite(1,*)'reduced'
1351
 
 
 
1350
                        
1352
1351
                        !from initial condition
1353
1352
                        call solve(pod_matrix%val, pod_rhs%val)
1354
1353
                        pod_coef_dt(:)=pod_rhs%val
1386
1385
                        print*,pod_coef
1387
1386
                     endif
1388
1387
 
1389
 
                  else !timestep.lg.1
 
1388
                  else !timestep.gt.1
1390
1389
                     ewrite(1,*)'timestep=',timestep
1391
 
 
 
1390
                     
1392
1391
                     POD_velocity=>extract_vector_field(POD_state(1,1,istate), "Velocity")
1393
1392
                     POD_pressure=>extract_scalar_field(POD_state(1,2,istate), "Pressure")
1394
 
                    ! call allocate(pod_matrix, POD_velocity, POD_pressure)
1395
 
                    ! call allocate(pod_rhs, POD_velocity, POD_pressure)
1396
1393
 
1397
1394
                     !pod_coef is from the previous timestep, argument?
1398
1395
                     open(40,file='pod_coef')
1399
1396
                     read(40,*)(pod_coef(i),i=1,(u%dim+1)*size(POD_state,1))
1400
1397
                     close(40)
1401
 
        ! if(.not.reduced_model) then
1402
1398
                     open(unit=20,file='pod_matrix_snapmean')
1403
1399
                     read(20,*)((pod_matrix_snapmean(i,j),j=1,(u%dim+1)*size(POD_state,1)),i=1,(u%dim+1)*size(POD_state,1))
1404
1400
                     read(20,*)(pod_rhs_snapmean(i),i=1,(u%dim+1)*size(POD_state,1))
1407
1403
                     pod_matrix%val=pod_matrix_snapmean(:,:)
1408
1404
                     pod_rhs%val=pod_rhs_snapmean(:)
1409
1405
 
1410
 
!goto 100
1411
1406
                     do k=1,size(POD_state,1)
1412
1407
                        do d=1, u%dim
1413
1408
 
1425
1420
                        pod_rhs%val=pod_rhs%val+pod_coef(k+u%dim*size(POD_state,1))*pod_rhs_perturbed(:)
1426
1421
                     enddo
1427
1422
 
1428
 
!100 continue
1429
 
                     open(51,file='pod_matrix_4000')
1430
 
                     write(51,*)(pod_matrix%val)
1431
 
                     close(51)
1432
 
                     open(52,file='pod_rhs_4000')
1433
 
                     write(52,*)(pod_rhs%val)
1434
 
                     close(52)
1435
 
               !  endif !.not.reduced_model
1436
1423
                       if(.not.reduced_model) then
1437
1424
                        pod_matrix%val=0.0
1438
1425
                        pod_matrix_mass%val=0.0
1466
1453
        
1467
1454
                     enddo  !k=1,size(POD_state,1)
1468
1455
                     
1469
 
                     open(51,file='pod_matrix_4000_2')
1470
 
                     write(51,*)(pod_matrix%val)
1471
 
                     close(51)
1472
 
                     open(52,file='pod_rhs_4000_2')
1473
 
                     write(52,*)(pod_rhs%val)
1474
 
                     close(52)
1475
 
                     !stop 45
1476
1456
                    endif
1477
1457
 
1478
 
                     if(petrov) then
1479
 
  print *,"before POD_u"
1480
 
    POD_u=>extract_vector_field(POD_state(1,1,istate), "Velocity")
1481
 
    print *,"before POD_u_scalar"
1482
 
    !POD_u_scalar=extract_scalar_field(POD_u,1)    !!need modify
1483
 
    POD_num=size(POD_state)
1484
 
    print *,"before allocate(X_val"
1485
 
    allocate(X_val(POD_u%dim,ele_loc(POD_u,ele)))
1486
 
    print *,"before allocate JJ"
1487
 
    allocate(JJ(POD_u%dim,mesh_dim(POD_u)))
1488
 
    !allocate(pod_coef(POD_u%dim*POD_num+POD_num))
1489
 
    pod_umesh => extract_mesh(POD_state(1,1,istate), "Mesh", stat)
1490
 
    TOTELE=pod_umesh%elements
1491
 
    NLOC=pod_umesh%shape%loc
1492
 
    NGI=pod_umesh%shape%ngi
1493
 
    print *,"before allocate rhs_old"
1494
 
   call allocate(pod_rhs_old, POD_u, POD_pressure)
1495
 
 
1496
 
    xf_shape=>ele_shape(POD_u,1)
1497
 
    x_shape=>ele_shape(POD_u,1)
1498
 
     print *, "timestep,, 2+++"
1499
 
do j=1, POD_num
1500
 
 !POD_u_scalar=>extract_scalar_field(POD_state(j), "PODPressure")
1501
 
       !POD_u_scalar=extract_scalar_field(POD_u,j)
1502
 
       X_val=ele_val(POD_u, j)
1503
 
       xf_shape=>ele_shape(Pod_u,j)  !!
1504
 
       x_shape=>ele_shape(POD_u,j)!!
1505
 
end do
1506
 
                        !!new added for petro-galerkin     
 
1458
                    !! Petrov-Galerkin POD
 
1459
                    if(petrov) then
 
1460
                       POD_u=>extract_vector_field(POD_state(1,1,istate), "Velocity")
 
1461
                       POD_num=size(POD_state)
 
1462
                       allocate(X_val(POD_u%dim,ele_loc(POD_u,ele)))
 
1463
                       allocate(JJ(POD_u%dim,mesh_dim(POD_u)))
 
1464
                       pod_umesh => extract_mesh(POD_state(1,1,istate), "Mesh", stat)
 
1465
                       TOTELE=pod_umesh%elements
 
1466
                       NLOC=pod_umesh%shape%loc
 
1467
                       NGI=pod_umesh%shape%ngi
 
1468
                       call allocate(pod_rhs_old, POD_u, POD_pressure)
 
1469
                       
 
1470
                       xf_shape=>ele_shape(POD_u,1)
 
1471
                       x_shape=>ele_shape(POD_u,1)
 
1472
                       do j=1, POD_num
 
1473
                          X_val=ele_val(POD_u, j)
 
1474
                          xf_shape=>ele_shape(Pod_u,j)  !!
 
1475
                          x_shape=>ele_shape(POD_u,j)!!
 
1476
                       end do
 
1477
                       !!new added for petro-galerkin     
1507
1478
       
1508
 
      
1509
 
       ! allocate(pod_sol_velocity(u_nodes,POD_velocity%dim))
1510
 
         allocate(theta_pet(NGI))
1511
 
         allocate(leftsvd(TOTELE*NLOC+NLOC,nsvd))
1512
 
         allocate(leftsvd_gi(NGI,nsvd))
1513
 
         allocate(leftsvd_x_gi(NGI,nsvd))
1514
 
         allocate(smean_gi(NGI))
1515
 
         allocate(AMAT_pod(nsvd,nsvd))
1516
 
         allocate(KMAT_pod(nsvd,nsvd))
1517
 
         allocate(KB_pod(nsvd)) 
1518
 
         allocate(b_pod(nsvd)) 
1519
 
         allocate(PSI(TOTELE*NLOC+NLOC))
1520
 
         allocate(PSI_GI(NGI)) 
1521
 
         allocate(PSI_OLD_GI(NGI)) 
1522
 
         allocate(psi_old(TOTELE*NLOC+NLOC))
1523
 
         allocate(GRADX(NGI))
1524
 
         allocate(GRADT(NGI))
1525
 
         allocate(DIFFGI_pod(NGI))
1526
 
         allocate(detwei(NGI))
1527
 
         allocate(res(NGI))
1528
 
       !  noloc=x_shape%loc
1529
 
       !  nogas=x_shape%ngi
1530
 
       !  allocate(dshape(noloc,nogas,POD_u%dim))
1531
 
         allocate(dshape(x_shape%loc,x_shape%ngi,POD_u%dim))  !number of nodes, NO. of gussian
 
1479
                       
 
1480
                       ! allocate(pod_sol_velocity(u_nodes,POD_velocity%dim))
 
1481
                       allocate(theta_pet(NGI))
 
1482
                       allocate(leftsvd(TOTELE*NLOC+NLOC,nsvd))
 
1483
                       allocate(leftsvd_gi(NGI,nsvd))
 
1484
                       allocate(leftsvd_x_gi(NGI,nsvd))
 
1485
                       allocate(smean_gi(NGI))
 
1486
                       allocate(AMAT_pod(nsvd,nsvd))
 
1487
                       allocate(KMAT_pod(nsvd,nsvd))
 
1488
                       allocate(KB_pod(nsvd)) 
 
1489
                       allocate(b_pod(nsvd)) 
 
1490
                       allocate(PSI(TOTELE*NLOC+NLOC))
 
1491
                       allocate(PSI_GI(NGI)) 
 
1492
                       allocate(PSI_OLD_GI(NGI)) 
 
1493
                       allocate(psi_old(TOTELE*NLOC+NLOC))
 
1494
                       allocate(GRADX(NGI))
 
1495
                       allocate(GRADT(NGI))
 
1496
                       allocate(DIFFGI_pod(NGI))
 
1497
                       allocate(detwei(NGI))
 
1498
                       allocate(res(NGI))
 
1499
                       allocate(dshape(x_shape%loc,x_shape%ngi,POD_u%dim))  !number of nodes, NO. of gussian
1532
1500
        
1533
1501
 
1534
 
     ! new added for petro-galerkin Residual, see  Non-Linear Petrov-Galerkin Methods for
1535
 
     ! Hyperbolic Equations and Discontinuous Finite Element Methods .. eqn 102
1536
 
     ! get the result of Residual value"
1537
 
     
1538
 
     pod_rhs_old%val=pod_rhs%val 
1539
 
     call solve(pod_matrix%val, pod_rhs%val)
1540
 
    ! pod_coef=pod_rhs%val
1541
 
     pod_rhs%val=pod_rhs%val-pod_rhs_old%val
1542
 
     call SMLINN(pod_matrix%val,RES,pod_rhs%val,GI,GI) ! get the RES
1543
 
 
1544
 
           DO ELE=1,TOTELE
1545
 
           leftsvd_gi = 0.0
1546
 
           leftsvd_x_gi = 0.0
1547
 
           smean_gi =0.0
1548
 
 
1549
 
          call transform_to_physical(POD_u, ele,  ele_shape(POD_u,1), detwei=detwei, dshape=dshape)
1550
 
           DO GI=1,NGI        ! guassian points' number
1551
 
             DO JLOC=1,NLOC   ! NLOC points in one element
1552
 
               GLOBJ=(ELE-1)*NLOC+JLOC
 
1502
                       ! new added for petro-galerkin Residual, see  Non-Linear Petrov-Galerkin Methods for
 
1503
                       ! Hyperbolic Equations and Discontinuous Finite Element Methods .. eqn 102
 
1504
                       ! get the result of Residual value"
 
1505
                       
 
1506
                       pod_rhs_old%val=pod_rhs%val 
 
1507
                       call solve(pod_matrix%val, pod_rhs%val)
 
1508
                       ! pod_coef=pod_rhs%val
 
1509
                       pod_rhs%val=pod_rhs%val-pod_rhs_old%val
 
1510
                       call SMLINN(pod_matrix%val,RES,pod_rhs%val,GI,GI) ! get the RES
 
1511
 
 
1512
                       DO ELE=1,TOTELE
 
1513
                          leftsvd_gi = 0.0
 
1514
                          leftsvd_x_gi = 0.0
 
1515
                          smean_gi =0.0
 
1516
                          
 
1517
                          call transform_to_physical(POD_u, ele,  ele_shape(POD_u,1), detwei=detwei, dshape=dshape)
 
1518
                          DO GI=1,NGI        ! guassian points' number
 
1519
                             DO JLOC=1,NLOC   ! NLOC points in one element
 
1520
                                GLOBJ=(ELE-1)*NLOC+JLOC
 
1521
                                
 
1522
                                do isvd =1,nsvd
 
1523
                                   leftsvd(GLOBJ,isvd)=POD_u%val(GLOBJ,isvd)
 
1524
                                   leftsvd_gi(GI,isvd)=leftsvd_gi(GI,isvd)+xf_shape%dn(isvd,GI,1)*leftsvd(GLOBJ,isvd) 
 
1525
                                   !N shape function, leftsvd: podbasis, N(JLOC,GI): ele_shape(nu, ele)
 
1526
                                   leftsvd_x_gi(GI,isvd)=leftsvd_gi(GI,isvd)+dshape(JLOC,GI,1)*leftsvd(GLOBJ,isvd)
 
1527
                                enddo
 
1528
                                !smean_gi(gi)=smean_gi(gi)+ 1SMLINN(A,X,B,NMX,N)*smean(GLOBJ)  !need modification
 
1529
                                !smean_gi(gi)=smean_gi(gi)+ 1*smean(GLOBJ)
 
1530
                             END DO
 
1531
                          END DO
 
1532
                       END DO
 
1533
 
 
1534
                       PSI_GI=0.0   !  quadrature point
 
1535
                       PSI_OLD_GI=0.0
 
1536
                       
 
1537
                       DO GI=1,NGI
 
1538
                          DO JLOC=1,NLOC   
 
1539
                             GLOBJ=(ELE-1)*NLOC+JLOC
 
1540
                             PSI_GI(GI)=PSI_GI(GI)+xf_shape%dn(isvd,GI,1)*PSI(GLOBJ)
 
1541
                             PSI_OLD_GI(GI)=PSI_OLD_GI(GI)+xf_shape%dn(isvd,GI,1)*PSI_OLD(GLOBJ)
 
1542
                          END DO
 
1543
                       END DO
 
1544
                       
 
1545
                       
 
1546
                       ! calculate the GRADX and GRADT    quadrature point
 
1547
                       GRADX=0.0
 
1548
                       GRADT=0.0
 
1549
                       DO GI=1,NGI      !   ?? NGI and NLOC need input 
 
1550
                          DO JLOC=1,NLOC
 
1551
                             GLOBJ=(ELE-1)*NLOC+JLOC
 
1552
                             PSI_theta=theta_pet(GI)*PSI(GLOBJ)+(1.-theta_pet(GI))*PSI_OLD(GLOBJ)
 
1553
                             
 
1554
                             GRADX(GI)=GRADX(GI)+dshape(JLOC,GI,1)*PSI_theta
 
1555
                             GRADT(GI)=GRADT(GI)+xf_shape%dn(isvd,GI,1)*(PSI(GLOBJ)-PSI_OLD(GLOBJ))/DT1
 
1556
                          END DO
 
1557
                       END DO
 
1558
 
 
1559
                       !  get the result of  AX_STAR cos P_STAR_POD need it 
 
1560
                       DO GI=1,NGI
 
1561
                          
 
1562
                          A_STAR_COEF=(GRADX(GI)+GRADT(GI)*1.) &
 
1563
                               /MAX(1.E-6,GRADX(GI)**2+1.*GRADT(GI)**2)
 
1564
                          !  AT_STAR=A_STAR_COEF*GRADT(GI)
 
1565
                          AX_STAR=A_STAR_COEF*GRADX(GI)
 
1566
                       enddo
 
1567
                       do gi =1,ngi
 
1568
                          do isvd = 1,nsvd
 
1569
                             if(isvd.eq.1) then
 
1570
                                P_STAR_POD = 0.25/abs(AX_STAR*leftsvd_x_gi(GI,isvd))
 
1571
                             else
 
1572
                                P_STAR_POD = 0.25/max(P_STAR_POD,abs(AX_STAR*leftsvd_x_gi(GI,isvd)))
 
1573
                             endif
 
1574
                          enddo
 
1575
                          DIFFGI_pod(GI) = 1.0e-4*P_STAR_POD*RES(GI)**2   &
 
1576
                               / MAX(1.E-6,GRADX(GI)**2+GRADT(GI)**2)!!v_pod (eqn.100)
 
1577
                       enddo
 
1578
                       
 
1579
                       do isvd =1,nsvd
 
1580
                          do jsvd=1,nsvd
 
1581
                             do gi =1,ngi
 
1582
                                KMAT_pod(isvd,jsvd) = KMAT_pod(isvd,jsvd) &
 
1583
                                     + detwei(GI)*leftsvd_x_gi(GI,isvd)*DIFFGI_pod(GI)*leftsvd_x_gi(GI,jsvd)
 
1584
                                ! (eqn.99 left side,leftsvd_x_gi:delta N)
 
1585
                             enddo
 
1586
                          enddo
 
1587
                       enddo
 
1588
                       
 
1589
 
 
1590
                       !!call solver
 
1591
                       
 
1592
                       
 
1593
                       
 
1594
                       !pod_rhs%val=0.001*pod_rhs%val
 
1595
                       ! add K to (MT AM + K) fi *pod= MT b. 
 
1596
                       
 
1597
                       do isvd =1,nsvd
 
1598
                          do jsvd=1,nsvd
 
1599
                             pod_matrix%val(isvd,jsvd) = pod_matrix%val(isvd,jsvd)+KMAT_pod(isvd,jsvd)
 
1600
                          enddo
 
1601
                          !b_pod(isvd) = b_pod(isvd)-KB_pod(isvd)
 
1602
                       enddo
 
1603
  
 
1604
                       call solve(pod_matrix%val,pod_rhs_old%val) 
 
1605
                       pod_rhs%val=pod_rhs_old%val
 
1606
                       deallocate(theta_pet)
 
1607
                       deallocate(leftsvd)
 
1608
                       deallocate(leftsvd_gi)
 
1609
                       deallocate(leftsvd_x_gi)
 
1610
                       deallocate(smean_gi)
 
1611
                       deallocate(AMAT_pod)
 
1612
                       deallocate(KMAT_pod)
 
1613
                       deallocate(KB_pod) 
 
1614
                       deallocate(b_pod) 
 
1615
                       deallocate(PSI)
 
1616
                       deallocate(PSI_GI) 
 
1617
                       deallocate(PSI_OLD_GI) 
 
1618
                       deallocate(psi_old)
 
1619
                       deallocate(GRADX)
 
1620
                       deallocate(GRADT)
 
1621
                       deallocate(DIFFGI_pod)
 
1622
                       deallocate(detwei)
 
1623
                       deallocate(res)
 
1624
                       deallocate(dshape) 
 
1625
                    else 
 
1626
                       call solve(pod_matrix%val,pod_rhs%val)
 
1627
                       
 
1628
                    endif ! end petrov 
 
1629
 
 
1630
                    ! call solve(pod_matrix%val, pod_rhs%val)
 
1631
                    pod_coef_dt=0.0
 
1632
                    pod_coef_dt=pod_rhs%val
 
1633
                    pod_rhs%val=0.0
 
1634
                    
 
1635
                    pod_coef(:)=pod_coef(:)+dt*pod_coef_dt(:)
 
1636
                    !pod_coef(:)=pod_coef_dt(:)
 
1637
                    
 
1638
                    !save pod_coef, rewrite every timestep
 
1639
                    open(40,file='pod_coef')
 
1640
                    write(40,*)(pod_coef(i),i=1,(u%dim+1)*size(POD_state,1))
 
1641
                    close(40)
 
1642
                    print*,pod_coef
 
1643
 
 
1644
                    call project_full(delta_u, delta_p, pod_sol_velocity, pod_sol_pressure, POD_state(:,:,istate), pod_coef)
 
1645
                    
 
1646
                    
 
1647
                 endif ! timestep          
 
1648
                 
 
1649
                 !1000 continue
 
1650
                 
 
1651
                 snapmean_velocity=>extract_vector_field(POD_state(1,1,istate),"SnapmeanVelocity")
 
1652
                 snapmean_pressure=>extract_Scalar_field(POD_state(1,2,istate),"SnapmeanPressure")
 
1653
                 !call addto(u, delta_u, dt)
 
1654
                 u%val=snapmean_velocity%val
 
1655
                 call addto(u, delta_u)
 
1656
                 p%val=snapmean_pressure%val
 
1657
                 call addto(p, delta_p)
 
1658
                 
 
1659
                 call deallocate(delta_p)
 
1660
                 call deallocate(delta_u)
 
1661
                 
 
1662
                 deallocate(pod_matrix_snapmean)
 
1663
                 deallocate(pod_rhs_snapmean)
 
1664
                 deallocate(pod_matrix_perturbed)
 
1665
                 deallocate(pod_matrix_adv_perturbed)
 
1666
                 deallocate(pod_rhs_perturbed)
 
1667
                 deallocate(pod_rhs_adv_perturbed)
 
1668
                 deallocate(pod_ct_m)
 
1669
                 call deallocate(pod_matrix)
 
1670
                 call deallocate(pod_matrix_mass)
 
1671
                 call deallocate(pod_rhs)
 
1672
                 
 
1673
                 
 
1674
              endif ! prognostic velocity
1553
1675
              
1554
 
               do isvd =1,nsvd
1555
 
                  leftsvd(GLOBJ,isvd)=POD_u%val(GLOBJ,isvd)
1556
 
                  leftsvd_gi(GI,isvd)=leftsvd_gi(GI,isvd)+xf_shape%dn(isvd,GI,1)*leftsvd(GLOBJ,isvd) 
1557
 
                      !N shape function, leftsvd: podbasis, N(JLOC,GI): ele_shape(nu, ele)
1558
 
                  leftsvd_x_gi(GI,isvd)=leftsvd_gi(GI,isvd)+dshape(JLOC,GI,1)*leftsvd(GLOBJ,isvd)
1559
 
               enddo
1560
 
                 !smean_gi(gi)=smean_gi(gi)+ 1SMLINN(A,X,B,NMX,N)*smean(GLOBJ)  !need modification
1561
 
                 !smean_gi(gi)=smean_gi(gi)+ 1*smean(GLOBJ)
1562
 
             END DO
1563
 
           END DO
1564
 
           END DO
1565
 
 
1566
 
           PSI_GI=0.0   !  quadrature point
1567
 
           PSI_OLD_GI=0.0
1568
 
           
1569
 
           DO GI=1,NGI
1570
 
             DO JLOC=1,NLOC   
1571
 
               GLOBJ=(ELE-1)*NLOC+JLOC
1572
 
               PSI_GI(GI)=PSI_GI(GI)+xf_shape%dn(isvd,GI,1)*PSI(GLOBJ)
1573
 
               PSI_OLD_GI(GI)=PSI_OLD_GI(GI)+xf_shape%dn(isvd,GI,1)*PSI_OLD(GLOBJ)
1574
 
             END DO
1575
 
           END DO
1576
 
        
1577
 
 
1578
 
 ! calculate the GRADX and GRADT    quadrature point
1579
 
           GRADX=0.0
1580
 
           GRADT=0.0
1581
 
           DO GI=1,NGI      !   ?? NGI and NLOC need input 
1582
 
             DO JLOC=1,NLOC
1583
 
               GLOBJ=(ELE-1)*NLOC+JLOC
1584
 
               PSI_theta=theta_pet(GI)*PSI(GLOBJ)+(1.-theta_pet(GI))*PSI_OLD(GLOBJ)
1585
 
            
1586
 
               GRADX(GI)=GRADX(GI)+dshape(JLOC,GI,1)*PSI_theta
1587
 
               GRADT(GI)=GRADT(GI)+xf_shape%dn(isvd,GI,1)*(PSI(GLOBJ)-PSI_OLD(GLOBJ))/DT1
1588
 
             END DO
1589
 
           END DO
1590
 
 
1591
 
!  get the result of  AX_STAR cos P_STAR_POD need it 
1592
 
         DO GI=1,NGI
1593
 
 
1594
 
             A_STAR_COEF=(GRADX(GI)+GRADT(GI)*1.) &
1595
 
                         /MAX(1.E-6,GRADX(GI)**2+1.*GRADT(GI)**2)
1596
 
            !  AT_STAR=A_STAR_COEF*GRADT(GI)
1597
 
             AX_STAR=A_STAR_COEF*GRADX(GI)
1598
 
           enddo
1599
 
         do gi =1,ngi
1600
 
             do isvd = 1,nsvd
1601
 
                if(isvd.eq.1) then
1602
 
                   P_STAR_POD = 0.25/abs(AX_STAR*leftsvd_x_gi(GI,isvd))
1603
 
                else
1604
 
                   P_STAR_POD = 0.25/max(P_STAR_POD,abs(AX_STAR*leftsvd_x_gi(GI,isvd)))
1605
 
                endif
1606
 
             enddo
1607
 
             DIFFGI_pod(GI) = 1.0e-4*P_STAR_POD*RES(GI)**2   &
1608
 
                            / MAX(1.E-6,GRADX(GI)**2+GRADT(GI)**2)!!v_pod (eqn.100)
1609
 
          enddo
1610
 
 
1611
 
          do isvd =1,nsvd
1612
 
             do jsvd=1,nsvd
1613
 
                do gi =1,ngi
1614
 
                   KMAT_pod(isvd,jsvd) = KMAT_pod(isvd,jsvd) &
1615
 
                                       + detwei(GI)*leftsvd_x_gi(GI,isvd)*DIFFGI_pod(GI)*leftsvd_x_gi(GI,jsvd)
1616
 
                    ! (eqn.99 left side,leftsvd_x_gi:delta N)
1617
 
                enddo
1618
 
             enddo
1619
 
          enddo 
1620
 
 
1621
 
 
1622
 
!!call solver
1623
 
 
1624
 
 
1625
 
   
1626
 
     !pod_rhs%val=0.001*pod_rhs%val
1627
 
     ! add K to (MT AM + K) fi *pod= MT b. 
1628
 
 
1629
 
       do isvd =1,nsvd
1630
 
               do jsvd=1,nsvd
1631
 
                  pod_matrix%val(isvd,jsvd) = pod_matrix%val(isvd,jsvd)+KMAT_pod(isvd,jsvd)
1632
 
               enddo
1633
 
                  !b_pod(isvd) = b_pod(isvd)-KB_pod(isvd)
1634
 
            enddo
1635
 
  
1636
 
       call solve(pod_matrix%val,pod_rhs_old%val) 
1637
 
       pod_rhs%val=pod_rhs_old%val
1638
 
          deallocate(theta_pet)
1639
 
          deallocate(leftsvd)
1640
 
          deallocate(leftsvd_gi)
1641
 
          deallocate(leftsvd_x_gi)
1642
 
          deallocate(smean_gi)
1643
 
          deallocate(AMAT_pod)
1644
 
          deallocate(KMAT_pod)
1645
 
          deallocate(KB_pod) 
1646
 
          deallocate(b_pod) 
1647
 
          deallocate(PSI)
1648
 
          deallocate(PSI_GI) 
1649
 
          deallocate(PSI_OLD_GI) 
1650
 
          deallocate(psi_old)
1651
 
          deallocate(GRADX)
1652
 
          deallocate(GRADT)
1653
 
          deallocate(DIFFGI_pod)
1654
 
          deallocate(detwei)
1655
 
          deallocate(res)
1656
 
          deallocate(dshape) 
1657
 
  else 
1658
 
       call solve(pod_matrix%val,pod_rhs%val)
1659
 
       
1660
 
 endif ! end petrov 
1661
 
                    ! call solve(pod_matrix%val, pod_rhs%val)
1662
 
                     pod_coef_dt=0.0
1663
 
                     pod_coef_dt=pod_rhs%val
1664
 
                     pod_rhs%val=0.0
1665
 
 
1666
 
                     pod_coef(:)=pod_coef(:)+dt*pod_coef_dt(:)
1667
 
                     !pod_coef(:)=pod_coef_dt(:)
1668
 
 
1669
 
                     !save pod_coef, rewrite every timestep
1670
 
                     open(40,file='pod_coef')
1671
 
                     write(40,*)(pod_coef(i),i=1,(u%dim+1)*size(POD_state,1))
1672
 
                     close(40)
1673
 
                     print*,pod_coef
1674
 
 
1675
 
                     call project_full(delta_u, delta_p, pod_sol_velocity, pod_sol_pressure, POD_state(:,:,istate), pod_coef)
1676
 
 
1677
 
 
1678
 
                  endif ! timestep          
1679
 
 
1680
 
!1000 continue
1681
 
 
1682
 
                  snapmean_velocity=>extract_vector_field(POD_state(1,1,istate),"SnapmeanVelocity")
1683
 
                  snapmean_pressure=>extract_Scalar_field(POD_state(1,2,istate),"SnapmeanPressure")
1684
 
                  !call addto(u, delta_u, dt)
1685
 
                  u%val=snapmean_velocity%val
1686
 
                  call addto(u, delta_u)
1687
 
                  p%val=snapmean_pressure%val
1688
 
                  call addto(p, delta_p)
1689
 
 
1690
 
                  call deallocate(delta_p)
1691
 
                  call deallocate(delta_u)
1692
 
 
1693
 
                  deallocate(pod_matrix_snapmean)
1694
 
                  deallocate(pod_rhs_snapmean)
1695
 
                  deallocate(pod_matrix_perturbed)
1696
 
                  deallocate(pod_matrix_adv_perturbed)
1697
 
                  deallocate(pod_rhs_perturbed)
1698
 
                  deallocate(pod_rhs_adv_perturbed)
1699
 
                  deallocate(pod_ct_m)
1700
 
                  call deallocate(pod_matrix)
1701
 
                  call deallocate(pod_matrix_mass)
1702
 
                  call deallocate(pod_rhs)
1703
 
 
1704
 
 
1705
 
               endif ! prognostic velocity
1706
 
 
1707
 
            enddo reduced_model_loop
1708
 
            call profiler_toc("reduced_model_loop")
1709
 
 
1710
 
          !  endif ! its=nonlinear_iterations
1711
 
         end if ! end of 'if .not.reduced_model'
 
1676
           enddo reduced_model_loop
 
1677
           call profiler_toc("reduced_model_loop")
 
1678
 
 
1679
           !  endif ! its=nonlinear_iterations
 
1680
        end if ! end of 'if .not.reduced_model'
1712
1681
 
1713
1682
 
1714
1683
        if(.not.reduced_model .or. (reduced_model .and. timestep_check))then