~reducedmodelling/fluidity/ReducedModel

« back to all changes in this revision

Viewing changes to tools/IGW.F90

  • Committer: skramer
  • Date: 2010-11-16 15:37:38 UTC
  • Revision ID: svn-v4:5bf5533e-7014-46e3-b1bb-cce4b9d03719:trunk:2477
Changing the storage of vector_field from vfield%val(dim)%ptr(node) to vfield%val(dim,node).

This was a long desired and discussed change that should improve cache performance and will allow vector fields of more than 3 dimension (internally no vtk support). Just realised it could be done with a few simple seds.

Most important changes need:
- removal of wrap_field() for vector_fields.
- this required some surgery in vtk_write_fields, a horrible piece of code in any case
- traffic and mesh_movement were using wrap_field and now have to do a few more copies.
- halo_update() requires a contiguous buffer in memory, so a copy is now done for halo_update(vfield)
- the wrapped python vector fields available in python_state diagnostics had to be rewired to account for the new storage. This actually simplifies the wrapping a lot.

I realise this is a big commit which might give some trouble when merging in big changes (from for example a branch). The following little script I used might be handy for this. You might only want to apply it to files you have modified. Make sure you back up your entire local checkout before applying any of these:

#/bin/sh
FILES=*/*.F90 */tests/*.F90
# replace any %val(1)%ptr(2) where 2 contains a pair of ()
sed -i 's!%val(\([A-Za-z0-9%_\+\* -]*\))%ptr(\([^)]*([^)]*)[^)]*\))!%val(\1,\2)!g' ${FILES}
# replace any %val(1)%ptr(2) where 2 contains no closing )
sed -i 's!%val(\([A-Za-z0-9%_\+\* -]*\))%ptr(\([^)]*\))!%val(\1,\2)!g' ${FILES}
# replace any %val(1)%ptr not followed by (
sed -i 's!%val(\([A-Za-z0-9%_\+\* -]*\))%ptr!%val(\1,:)!g' ${FILES}
# now svn revert femtools/Sparse_Tools.F90
# as it gives false positives for dcsr_matrices
# luckily it doesn't use vector_fields, so we can just revert


Show diffs side-by-side

added added

removed removed

Lines of Context:
378
378
          max_its=1000, atol=1.0e-30, rtol=1.0e-14, &
379
379
          start_from_zero=.true.)
380
380
     tmpu = Ro*tmpu
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)
385
385
     tmpu = -Ro*tmpu
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")
388
388
  else
389
389
     inquire(file=trim(u_input),exist=file_exists)  
549
549
     ewrite(2,*) 'computing div u'
550
550
     
551
551
     ctu%val = 0.0
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
556
556
 
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
586
586
 
587
587
     ewrite(2,*) 'computing div u'
588
588
     
589
589
     ctu%val = 0.0
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
594
594
 
595
595
     ewrite(2,*) 'max(abs(ctu))=',maxval(ctu%val), minval(ctu%val)
601
601
     ewrite(2,*) 'computing curl of u'
602
602
 
603
603
     ctu%val = 0.0
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
608
608
 
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
632
632
 
633
633
     call vtk_write_fields(trim('projection'), &
634
634
          index=0, position=positions, &
640
640
  end if
641
641
 
642
642
  ctu%val = 0.0
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
647
647
  
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/), &
673
673
     rhs = 0.
674
674
 
675
675
     !u mass parts
676
 
     call mult(tmpu,u_mass_static,u%val(1)%ptr)
 
676
     call mult(tmpu,u_mass_static,u%val(1,:))
677
677
     !d/dt
678
678
     rhs(1:node_count(u))=rhs(1:node_count(u)) + tmpu
679
679
     !coriolis
681
681
          rhs(1+node_count(u):2*node_count(u))&
682
682
          + 0.5*dt*tmpu*Fr*Fr/Ro
683
683
 
684
 
     call mult(tmpu,u_mass_static,u%val(2)%ptr)
 
684
     call mult(tmpu,u_mass_static,u%val(2,:))
685
685
     !d/dt
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)) - &
700
700
             0.5*dt*tmpu
701
701
        
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)) &
705
705
             +0.5*dt*tmph
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)) &
709
709
             +0.5*dt*tmph
731
731
        ewrite(1,*) 'dumping t = ',t
732
732
 
733
733
        ctu%val = 0.0
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
738
738
 
739
739
        !call petsc_solve(ctu%val,h_mass_static, ctu%val, &
741
741
        !     iterations=number_of_iterations)
742
742
        
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/), &