Ignore:
Timestamp:
Mar 30, 2013, 2:44:52 PM (13 years ago)
Author:
Julian Iseringhausen <isering@…>
Children:
8180d8
Parents:
d13e27
git-author:
Julian Iseringhausen <isering@…> (06/11/12 14:02:16)
git-committer:
Julian Iseringhausen <isering@…> (03/30/13 14:44:52)
Message:

Open boundary conditions.

Conflicts:

lib/vmg/src/Makefile.am
lib/vmg/src/base/factory.cpp
lib/vmg/test/unit_test/library/dirichlet_fas_lr_mpi.cpp

File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/samples/discretization_poisson_fv.cpp

    rd13e27 rf57182  
    3737using namespace VMG;
    3838
     39void DiscretizationPoissonFV::InitDiscretizationPoissonFV()
     40{
     41  switch (order)
     42    {
     43    case 2:
     44      stencil.SetDiag(6.0);
     45      stencil.push_back(-1,  0,  0, -1.0);
     46      stencil.push_back( 1,  0,  0, -1.0);
     47      stencil.push_back( 0, -1,  0, -1.0);
     48      stencil.push_back( 0,  1,  0, -1.0);
     49      stencil.push_back( 0,  0, -1, -1.0);
     50      stencil.push_back( 0,  0,  1, -1.0);
     51      break;
     52    case 4:
     53      stencil.SetDiag(24.0/6.0);
     54      stencil.push_back(-1,  0,  0, -2.0/6.0);
     55      stencil.push_back( 1,  0,  0, -2.0/6.0);
     56      stencil.push_back( 0, -1,  0, -2.0/6.0);
     57      stencil.push_back( 0,  1,  0, -2.0/6.0);
     58      stencil.push_back( 0,  0, -1, -2.0/6.0);
     59      stencil.push_back( 0,  0,  1, -2.0/6.0);
     60      stencil.push_back(-1, -1,  0, -1.0/6.0);
     61      stencil.push_back(-1,  1,  0, -1.0/6.0);
     62      stencil.push_back( 1, -1,  0, -1.0/6.0);
     63      stencil.push_back( 1,  1,  0, -1.0/6.0);
     64      stencil.push_back(-1,  0, -1, -1.0/6.0);
     65      stencil.push_back(-1,  0,  1, -1.0/6.0);
     66      stencil.push_back( 1,  0, -1, -1.0/6.0);
     67      stencil.push_back( 1,  0,  1, -1.0/6.0);
     68      stencil.push_back( 0, -1, -1, -1.0/6.0);
     69      stencil.push_back( 0, -1,  1, -1.0/6.0);
     70      stencil.push_back( 0,  1, -1, -1.0/6.0);
     71      stencil.push_back( 0,  1,  1, -1.0/6.0);
     72      break;
     73    default:
     74      assert(0 != "vmg choose discretization order 2 or 4");
     75      break;
     76    }
     77}
     78
     79void DiscretizationPoissonFV::ModifyRightHandSide()
     80{
     81  if (order == 4) {
     82
     83    Grid& rhs = MG::GetRhsMaxLevel();
     84
     85    Stencil stencil(6.0/12.0);
     86    stencil.push_back(-1,  0,  0, 1.0/12.0);
     87    stencil.push_back( 1,  0,  0, 1.0/12.0);
     88    stencil.push_back( 0, -1,  0, 1.0/12.0);
     89    stencil.push_back( 0,  1,  0, 1.0/12.0);
     90    stencil.push_back( 0,  0, -1, 1.0/12.0);
     91    stencil.push_back( 0,  0,  1, 1.0/12.0);
     92
     93    stencil.Apply(rhs);
     94
     95  }
     96}
     97
    3998void DiscretizationPoissonFV::SetInnerBoundaryCompute(Grid& sol_f, Grid& rhs_f, Grid& sol_c) const
    4099{
    … …  
    48107  const Index b2_f = rhs_f.Local().SizeTotal() - 1;
    49108
    50   const Index begin_f = sol_f.Local().Begin();
    51   const Index end_f = sol_f.Local().End();
    52   const Index begin_c = sol_c.Local().FinerBegin();
     109  const Index& begin_f = sol_f.Local().Begin();
     110  const Index& end_f = sol_f.Local().End();
     111  const Index& begin_c = sol_c.Local().FinerBegin();
    53112
    54113  const vmg_float c_1_3 = 1.0 / 3.0;
    … …  
    56115  const vmg_float c_4_3 = 4.0 / 3.0;
    57116
     117  Helper::AssertVectorsEqual(sol_c.GetSpatialPos(sol_c.Local().FinerBegin()), sol_f.GetSpatialPos(sol_f.Local().Begin()));
     118  Helper::AssertVectorsEqual(sol_c.GetSpatialPos(sol_c.Local().FinerEnd()-1), sol_f.GetSpatialPos(sol_f.Local().End()-1));
     119
    58120  //
    59121  // X-direction
    60122  //
     123  i_f = begin_f; i_c = begin_c;
    61124  for (i_f.Y()=begin_f.Y(), i_c.Y()=begin_c.Y(); i_f.Y()<end_f.Y(); i_f.Y()+=2, ++i_c.Y())
    62125    for (i_f.Z()=begin_f.Z(), i_c.Z()=begin_c.Z(); i_f.Z()<end_f.Z(); i_f.Z()+=2, ++i_c.Z()) {
     126      Helper::AssertVectorsEqual(sol_c.GetSpatialPos(i_c), sol_f.GetSpatialPos(i_f));
    63127      rhs_f(b1_f.X(),i_f.Y(),i_f.Z()) = (sol_c(b1_c.X()-1,i_c.Y(),i_c.Z()) - sol_f(b1_f.X()+1,i_f.Y(),i_f.Z())) * h2_inv.X();
    64128      rhs_f(b2_f.X(),i_f.Y(),i_f.Z()) = (sol_c(b2_c.X()+1,i_c.Y(),i_c.Z()) - sol_f(b2_f.X()-1,i_f.Y(),i_f.Z())) * h2_inv.X();
    … …  
    106170  // Y-direction
    107171  //
    108 
     172  i_f = begin_f; i_c = begin_c;
    109173  for (i_f.X()=begin_f.X(), i_c.X()=begin_c.X(); i_f.X()<end_f.X(); i_f.X()+=2, ++i_c.X())
    110174    for (i_f.Z()=begin_f.Z(), i_c.Z()=begin_c.Z(); i_f.Z()<end_f.Z(); i_f.Z()+=2, ++i_c.Z()) {
     175      Helper::AssertVectorsEqual(sol_c.GetSpatialPos(i_c), sol_f.GetSpatialPos(i_f));
    111176      rhs_f(i_f.X(),b1_f.Y(),i_f.Z()) = (sol_c(i_c.X(),b1_c.Y()-1,i_c.Z()) - sol_f(i_f.X(),b1_f.Y()+1,i_f.Z())) * h2_inv.Y();
    112177      rhs_f(i_f.X(),b2_f.Y(),i_f.Z()) = (sol_c(i_c.X(),b2_c.Y()+1,i_c.Z()) - sol_f(i_f.X(),b2_f.Y()-1,i_f.Z())) * h2_inv.Y();
    … …  
    154219  // Z-direction
    155220  //
    156 
     221  i_f = begin_f; i_c = begin_c;
    157222  for (i_f.X()=begin_f.X(), i_c.X()=begin_c.X(); i_f.X()<end_f.X(); i_f.X()+=2, ++i_c.X())
    158223    for (i_f.Y()=begin_f.Y(), i_c.Y()=begin_c.Y(); i_f.Y()<end_f.Y(); i_f.Y()+=2, ++i_c.Y()) {
     224      Helper::AssertVectorsEqual(sol_c.GetSpatialPos(i_c), sol_f.GetSpatialPos(i_f));
    159225      rhs_f(i_f.X(),i_f.Y(),b1_f.Z()) = (sol_c(i_c.X(),i_c.Y(),b1_c.Z()-1) - sol_f(i_f.X(),i_f.Y(),b1_f.Z()+1)) * h2_inv.Z();
    160226      rhs_f(i_f.X(),i_f.Y(),b2_f.Z()) = (sol_c(i_c.X(),i_c.Y(),b2_c.Z()+1) - sol_f(i_f.X(),i_f.Y(),b2_f.Z()-1)) * h2_inv.Z();
Note: See TracChangeset for help on using the changeset viewer.