378
378
max_its=1000, atol=1.0e-30, rtol=1.0e-14, &
379
379
start_from_zero=.true.)
381
call petsc_solve(u%val(1)%ptr,u_mass_static,tmpu, &
381
call petsc_solve(u%val(1,:),u_mass_static,tmpu, &
382
382
option_path="/temporary/fix")
383
383
ewrite(2,*) 'getting balanced v'
384
384
call mult_T(tmpu,C1T_static,h%val)
386
call petsc_solve(u%val(2)%ptr,u_mass_static,tmpu, &
386
call petsc_solve(u%val(2,:),u_mass_static,tmpu, &
387
387
option_path="/temporary/fix")
389
389
inquire(file=trim(u_input),exist=file_exists)
549
549
ewrite(2,*) 'computing div u'
552
call mult(tmph,C1T_static,u%val(1)%ptr)
552
call mult(tmph,C1T_static,u%val(1,:))
553
553
ctu%val = ctu%val + tmph
554
call mult(tmph,C2T_static,u%val(2)%ptr)
554
call mult(tmph,C2T_static,u%val(2,:))
555
555
ctu%val = ctu%val + tmph
557
557
ewrite(2,*) 'max(abs(ctu))=',maxval(ctu%val), minval(ctu%val)
574
574
startfromzero=.true., checkconvergence=.true., &
575
575
iterations=number_of_iterations)
576
576
!call mult(tmpu,u_inverse_mass_static,tmpu)
577
u%val(1)%ptr = u%val(1)%ptr - tmpu
577
u%val(1,:) = u%val(1,:) - tmpu
578
578
ewrite(2,*) 'max u1 correction', maxval(abs(tmpu))
579
579
call mult_T(tmpu,C2T_static,tmph)
580
580
call petsc_solve(tmpu,u_mass_static,tmpu, &
582
582
iterations=number_of_iterations)
583
583
!call mult(tmpu,u_inverse_mass_static,tmpu)
584
584
ewrite(2,*) 'max u2 correction', maxval(abs(tmpu))
585
u%val(2)%ptr = u%val(2)%ptr - tmpu
585
u%val(2,:) = u%val(2,:) - tmpu
587
587
ewrite(2,*) 'computing div u'
590
call mult(tmph,C1T_static,u%val(1)%ptr)
590
call mult(tmph,C1T_static,u%val(1,:))
591
591
ctu%val = ctu%val + tmph
592
call mult(tmph,C2T_static,u%val(2)%ptr)
592
call mult(tmph,C2T_static,u%val(2,:))
593
593
ctu%val = ctu%val + tmph
595
595
ewrite(2,*) 'max(abs(ctu))=',maxval(ctu%val), minval(ctu%val)
601
601
ewrite(2,*) 'computing curl of u'
604
call mult(tmph,C2T_static,u%val(1)%ptr)
604
call mult(tmph,C2T_static,u%val(1,:))
605
605
ctu%val = ctu%val - tmph
606
call mult(tmph,C1T_static,u%val(2)%ptr)
606
call mult(tmph,C1T_static,u%val(2,:))
607
607
ctu%val = ctu%val + tmph
609
609
ewrite(2,*) 'max(abs(ctu))=',maxval(ctu%val), minval(ctu%val)
621
621
call petsc_solve(tmpu,u_mass_static,tmpu, &
622
622
startfromzero=.true., checkconvergence=.true., &
623
623
iterations=number_of_iterations)
624
ewrite(2,*) 'difference for u2', maxval(abs(tmpu-u%val(2)%ptr))
625
error_u%val(2)%ptr = tmpu
624
ewrite(2,*) 'difference for u2', maxval(abs(tmpu-u%val(2,:)))
625
error_u%val(2,:) = tmpu
626
626
call mult_T(tmpu,C2T_static,tmph)
627
627
call petsc_solve(tmpu,u_mass_static,tmpu, &
628
628
startfromzero=.true., checkconvergence=.true., &
629
629
iterations=number_of_iterations)
630
ewrite(2,*) 'difference for u1', maxval(abs(tmpu+u%val(1)%ptr))
631
error_u%val(1)%ptr = -tmpu
630
ewrite(2,*) 'difference for u1', maxval(abs(tmpu+u%val(1,:)))
631
error_u%val(1,:) = -tmpu
633
633
call vtk_write_fields(trim('projection'), &
634
634
index=0, position=positions, &
643
call mult(tmph,C1T_static,u%val(1)%ptr)
643
call mult(tmph,C1T_static,u%val(1,:))
644
644
ctu%val = ctu%val + tmph
645
call mult(tmph,C2T_static,u%val(2)%ptr)
645
call mult(tmph,C2T_static,u%val(2,:))
646
646
ctu%val = ctu%val + tmph
648
648
if(steady_state) then
649
initial_u%val(1)%ptr = u%val(1)%ptr
650
initial_u%val(2)%ptr = u%val(2)%ptr
651
error_u%val(1)%ptr = 0.
652
error_u%val(2)%ptr = 0.
649
initial_u%val(1,:) = u%val(1,:)
650
initial_u%val(2,:) = u%val(2,:)
651
error_u%val(1,:) = 0.
652
error_u%val(2,:) = 0.
653
653
call vtk_write_fields(trim(filename), &
654
654
index=0, position=positions, &
655
655
model=vtk_mesh, sfields=(/h,ctu/), &
676
call mult(tmpu,u_mass_static,u%val(1)%ptr)
676
call mult(tmpu,u_mass_static,u%val(1,:))
678
678
rhs(1:node_count(u))=rhs(1:node_count(u)) + tmpu
681
681
rhs(1+node_count(u):2*node_count(u))&
682
682
+ 0.5*dt*tmpu*Fr*Fr/Ro
684
call mult(tmpu,u_mass_static,u%val(2)%ptr)
684
call mult(tmpu,u_mass_static,u%val(2,:))
686
686
rhs(1+node_count(u):2*node_count(u))= &
687
687
rhs(1+node_count(u):2*node_count(u))&
699
699
rhs(1+node_count(u):2*node_count(u)) - &
702
call mult(tmph,C1T_static,u%val(1)%ptr)
702
call mult(tmph,C1T_static,u%val(1,:))
703
703
rhs(2*node_count(u)+1:2*node_count(u)+node_count(h)) = &
704
704
rhs(2*node_count(u)+1:2*node_count(u)+node_count(h)) &
706
call mult(tmph,C2T_static,u%val(2)%ptr)
706
call mult(tmph,C2T_static,u%val(2,:))
707
707
rhs(2*node_count(u)+1:2*node_count(u)+node_count(h)) = &
708
708
rhs(2*node_count(u)+1:2*node_count(u)+node_count(h)) &
731
731
ewrite(1,*) 'dumping t = ',t
734
call mult(tmph,C1T_static,u%val(1)%ptr)
734
call mult(tmph,C1T_static,u%val(1,:))
735
735
ctu%val = ctu%val + tmph
736
call mult(tmph,C2T_static,u%val(2)%ptr)
736
call mult(tmph,C2T_static,u%val(2,:))
737
737
ctu%val = ctu%val + tmph
739
739
!call petsc_solve(ctu%val,h_mass_static, ctu%val, &
741
741
! iterations=number_of_iterations)
743
743
if (steady_state) then
744
error_u%val(1)%ptr = u%val(1)%ptr - initial_u%val(1)%ptr
745
error_u%val(2)%ptr = u%val(2)%ptr - initial_u%val(2)%ptr
744
error_u%val(1,:) = u%val(1,:) - initial_u%val(1,:)
745
error_u%val(2,:) = u%val(2,:) - initial_u%val(2,:)
746
746
call vtk_write_fields(trim(filename), &
747
747
index=dump, position=positions, &
748
748
model=vtk_mesh, sfields=(/h,ctu/), &