| 60 | | and sometimes we may be implementing a function that is only defined for a certain region (like a clump object etc...) where the value for q outside of a given radius is unknown. In that case, we just want to add up changes to q for subcells within the region. So for example to calculate the x component of the potential along a given edge we would |
| 61 | | {{{ |
| 62 | | dx=levels(Info%level)%dx |
| 63 | | ddx=dx/N |
| 64 | | DO i=1,Info%mX(1)+1 |
| 65 | | pos(1)=Info%xBbounds(1,1)+(i-1)*dx |
| 66 | | DO j=1,Info%mx(2)+1 |
| 67 | | pos(2)=Info%xBbounds(1,1)+(i-1)*dx |
| 68 | | DO k=1,Info%mx(3)+1 |
| 69 | | pos(3)=Info%xBbounds(1,1)+(i-1)*dx |
| 70 | | drho=0d0 |
| 71 | | DO ii=1,N |
| 72 | | ppos(1)=pos(1)+(REAL(ii)-.5)*ddx |
| 73 | | DO jj=1,N |
| 74 | | ppos(2)=pos(2)+(REAL(jj)-.5)*ddx |
| 75 | | DO kk=1,N |
| 76 | | ppos(3)=pos(3)+(REAL(kk)-.5)*ddx |
| 77 | | drho=(rho_func(ppos)-Info%q(i,j,k,1)) |
| 78 | | END DO |
| 79 | | END DO |
| 80 | | END DO |
| 81 | | Info%q(i,j,k,1)=Info%q(i,j,k,1)+drho/N**nDim |
| 82 | | END DO |
| 83 | | END DO |
| 84 | | END DO |
| 85 | | }}} |
| 86 | | |
| 87 | | For face centered fields we could also subsample, but their is the divergence criterion to consider for B-fields. We could calculate the potential on a subgrid, and then take a bunch of curls on each face and add them up, but Stoke's theorem lets us just do the integral around the outside. So we really just need to subsample along the appropriate edge for each component of the vector potential. |
| | 60 | and sometimes we may be implementing a function that is only defined for a certain region (like within a certain radius etc...) where the value for q outside of a given radius is unknown and should be left unchanged. In that case, we just want to add up changes to q for subcells within the region. |
| | 89 | == Face centered quantities == |
| | 90 | * AstroBEAR stores magnetic fields as face averages in the array {{{Info%aux(i,j,k,:)}}} |
| | 91 | * Each cell is a face perpendicular to the component of the B-field with sides of length {{{dx=levels(Info%level)%dx}}} |
| | 92 | * What is more important is the location of the edges since those will store in temporary Vector Potentials which can then be used to update the aux fields. |
| | 93 | * The 'lower' end point of every component of the vector potential A(i,j,k,m) is at {{{Info%xbounds(:,1) + ((/i,j,k/)-1)*dx}}} |
| | 94 | * The 'upper' end point of each component of the vector potential A(i,j,k,m) is at {{{Info%xbounds(:,1) + ((/i,j,k/)+I(m,:))*dx}}} where {{{I(:,:)}}} is the identity matrix so it effectively adds 1 along the 'm' coordinate. |
| | 95 | |
| | 96 | === Sub-sampling === |
| | 97 | |
| | 98 | For face centered fields we could also subsample, but there is the divergence criterion to consider for B-fields. We could calculate the potential on a subgrid, and then take a bunch of curls on each face and add them up, but Stoke's theorem lets us just do the integral around the outside. So we really just need to subsample along the appropriate edge for each component of the vector potential. So for example to calculate the x component of the potential along a given edge we would |
| | 99 | {{{ |
| | 100 | dx=levels(Info%level)%dx |
| | 101 | ddx=dx/N |
| | 102 | DO i=1,Info%mX(1) |
| | 103 | pos(1)=Info%xBbounds(1,1)+(i-1)*dx |
| | 104 | DO j=1,Info%mx(2)+1 |
| | 105 | pos(2)=Info%xBbounds(1,1)+(i-1)*dx |
| | 106 | DO k=1,Info%mx(3)+1 |
| | 107 | pos(3)=Info%xBbounds(1,1)+(i-1)*dx |
| | 108 | A(i,j,k,1)=0d0 |
| | 109 | DO ii=1,N |
| | 110 | ppos(1)=pos(1)+(REAL(ii)-.5)*ddx |
| | 111 | A(i,j,k,1)=Ax(ppos) |
| | 112 | END DO |
| | 113 | A(i,j,k,1)=A(i,j,k,1)/N |
| | 114 | END DO |
| | 115 | END DO |
| | 116 | END DO |
| | 117 | }}} |
| | 118 | |
| | 119 | or we could do all three components of the potential like this |
| | 120 | {{{ |
| | 121 | dx=levels(Info%level)%dx |
| | 122 | ddx=dx/N |
| | 123 | DO m=1,3 |
| | 124 | p=1 |
| | 125 | p(m)=0 |
| | 126 | DO i=1,Info%mX(1)+p(1) |
| | 127 | pos(1)=Info%xBbounds(1,1)+(i-1)*dx |
| | 128 | DO j=1,Info%mx(2)+p(2) |
| | 129 | pos(2)=Info%xBbounds(1,1)+(i-1)*dx |
| | 130 | DO k=1,Info%mx(3)+p(3) |
| | 131 | pos(3)=Info%xBbounds(1,1)+(i-1)*dx |
| | 132 | A(i,j,k,m)=0d0 |
| | 133 | DO ii=1,N |
| | 134 | ppos(m)=pos(m)+(REAL(ii)-.5)*ddx |
| | 135 | temp(:)=VectorPotential(ppos) !This returns a 3-component vector for simplicity |
| | 136 | A(i,j,k,m)=A(i,j,k,m)+temp(m) !but we only retain the 'm' component. |
| | 137 | END DO |
| | 138 | A(i,j,k,m)=A(i,j,k,m)/N |
| | 139 | END DO |
| | 140 | END DO |
| | 141 | END DO |
| | 142 | END DO |
| | 143 | }}} |
| | 144 | |
| | 145 | There are many examples of loops like these in the various objects (ie clumps, outflows, disks, etc..) |
| | 146 | |
| | 147 | One could imagine writing functions like 'volume average', or 'edge average' that could replace those inner loops - that would require a user function to be passed. |
