Changeset 616 in Athena


Ignore:
Timestamp:
01/12/08 11:17:28 (6 years ago)
Author:
jstone
Message:

Fixed bug in V3.1; self-gravity not included in computation of cell-centered
emfs at n+1/2 in CTU+CT integrators.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • tags/version-3.1/src/integrate_3d_ctu.c

    r607 r616  
    14591459           - q2*(x2Flux[k  ][j+1][i  ].Mz - x2Flux[k][j][i].Mz) 
    14601460           - q3*(x3Flux[k+1][j  ][i  ].My - x3Flux[k][j][i].My); 
     1461 
     1462        M2 = pG->U[k][j][i].M2 
     1463           - q1*(x1Flux[k  ][j  ][i+1].My - x1Flux[k][j][i].My) 
     1464           - q2*(x2Flux[k  ][j+1][i  ].Mx - x2Flux[k][j][i].Mx) 
     1465           - q3*(x3Flux[k+1][j  ][i  ].Mz - x3Flux[k][j][i].Mz); 
     1466 
     1467        M3 = pG->U[k][j][i].M3 
     1468           - q1*(x1Flux[k  ][j  ][i+1].Mz - x1Flux[k][j][i].Mz) 
     1469           - q2*(x2Flux[k  ][j+1][i  ].My - x2Flux[k][j][i].My) 
     1470           - q3*(x3Flux[k+1][j  ][i  ].Mx - x3Flux[k][j][i].Mx); 
     1471 
     1472/* Add source terms for fixed gravitational potential */ 
    14611473        if (StaticGravPot != NULL){ 
    14621474          phir = (*StaticGravPot)((x1+0.5*pG->dx1),x2,x3); 
    14631475          phil = (*StaticGravPot)((x1-0.5*pG->dx1),x2,x3); 
    14641476          M1 -= q1*(phir-phil)*pG->U[k][j][i].d; 
    1465         } 
    1466  
    1467         M2 = pG->U[k][j][i].M2 
    1468            - q1*(x1Flux[k  ][j  ][i+1].My - x1Flux[k][j][i].My) 
    1469            - q2*(x2Flux[k  ][j+1][i  ].Mx - x2Flux[k][j][i].Mx) 
    1470            - q3*(x3Flux[k+1][j  ][i  ].Mz - x3Flux[k][j][i].Mz); 
    1471         if (StaticGravPot != NULL){ 
     1477 
    14721478          phir = (*StaticGravPot)(x1,(x2+0.5*pG->dx2),x3); 
    14731479          phil = (*StaticGravPot)(x1,(x2-0.5*pG->dx2),x3); 
    14741480          M2 -= q2*(phir-phil)*pG->U[k][j][i].d; 
    1475         } 
    1476  
    1477         M3 = pG->U[k][j][i].M3 
    1478            - q1*(x1Flux[k  ][j  ][i+1].Mz - x1Flux[k][j][i].Mz) 
    1479            - q2*(x2Flux[k  ][j+1][i  ].My - x2Flux[k][j][i].My) 
    1480            - q3*(x3Flux[k+1][j  ][i  ].Mx - x3Flux[k][j][i].Mx); 
    1481         if (StaticGravPot != NULL){ 
     1481 
    14821482          phir = (*StaticGravPot)(x1,x2,(x3+0.5*pG->dx3)); 
    14831483          phil = (*StaticGravPot)(x1,x2,(x3-0.5*pG->dx3)); 
    14841484          M3 -= q3*(phir-phil)*pG->U[k][j][i].d; 
    14851485        } 
     1486 
     1487/* Add source terms due to self-gravity  */ 
     1488#ifdef SELF_GRAVITY 
     1489        phir = 0.5*(pG->Phi[k][j][i] + pG->Phi[k][j][i+1]); 
     1490        phil = 0.5*(pG->Phi[k][j][i] + pG->Phi[k][j][i-1]); 
     1491        M1 -= q1*(phir-phil)*pG->U[k][j][i].d; 
     1492 
     1493        phir = 0.5*(pG->Phi[k][j][i] + pG->Phi[k][j+1][i]); 
     1494        phil = 0.5*(pG->Phi[k][j][i] + pG->Phi[k][j-1][i]); 
     1495        M2 -= q2*(phir-phil)*pG->U[k][j][i].d; 
     1496 
     1497        phir = 0.5*(pG->Phi[k][j][i] + pG->Phi[k+1][j][i]); 
     1498        phil = 0.5*(pG->Phi[k][j][i] + pG->Phi[k-1][j][i]); 
     1499        M3 -= q3*(phir-phil)*pG->U[k][j][i].d; 
     1500#endif /* SELF_GRAVITY */ 
    14861501 
    14871502        B1c = 0.5*(B1_x1Face[k][j][i] + B1_x1Face[k  ][j  ][i+1]); 
Note: See TracChangeset for help on using the changeset viewer.