/** * \file Solver.C * \brief Solver class definition. */ #include "Solver.h" #include "time.h" using namespace std; const real rIonosphere = RION*1.e-2; //!< in meters //Range all; /********************************************************************//** * \author Slava Merkin (vgm at bu.edu) * \since 06-2007 * * \param[in] LOW_LAT_BOUND: Colatitude of the solver low latitude * boundary in degrees. Requirements: * - 30 < LOW_LAT_BOUND < 60 * - Divisible by latitudinal grid step (LAT_GRID) * * \param[in] LAT_GRID: Latitudinal grid step in degrees. Requirements: * - LAT_GRID <= 5 degrees * - Divides evenly into low latitude boundary (LOW_LAT_BOUND) * - 180 % LAT_GRID == 0 and the result should be an even number * * \param[in] LON_GRID: Longitudinal Grid Step. Requirements: * - LON_GRID <= 5 degrees * - 360 % LON_GRID ==0 and the result should be an even number. * * The constructor first checks whether the parameters satisfy the * requirements listed in const.h. Then the necessary coordinate * arrays are dimensioned (nLat,nLon) and solverSphericalGrid is * created by calling generateSpherical2DGrid() ************************************************************************/ Solver::Solver(const real LOW_LAT_BOUND, const real LAT_GRID, const real LON_GRID) { /* if ( ( lowLatBound < 30 ) or ( lowLatBound > 60 ) or ( lowLatBound % latGrid != 0) or ( 360 % lonGrid != 0) or ( latGrid > 5 ) or ( lonGrid > 5 ) or ( (180 % latGrid) % 2 != 0) or ( (360 % lonGrid) % 2 != 0) ) { cout << "The solver grid parameters in const.h are incorrect. See const.h documentation."<< endl; exit(1);}*/ lowLatBound = LOW_LAT_BOUND; latGrid = LAT_GRID; lonGrid = LON_GRID; nLat = lowLatBound/latGrid+1; nLon = 360./lonGrid+1; positions.redim(nLat*nLon,2); x.redim(nLat,nLon); y.redim(nLat,nLon); latitude.redim(nLat,nLon); longitude.redim(nLat,nLon); realArray lon(nLon), lat(nLat); lon.seqAdd(0.,lonGrid); lat.seqAdd(0.,latGrid); lon = lon/180.*PI; lat = lat/180.*PI; // ===== The part below is the same for all constructors ========== Range NLAT(0,nLat-1); for (int k = 0; k < nLon; k++) { latitude(all,k) = lat; longitude(all,k) = lon(k); x(all,k)=sin(lat)*cos(lon(k)); y(all,k)=sin(lat)*sin(lon(k)); positions(NLAT+k*nLat,0) = x(all,k); positions(NLAT+k*nLat,1) = y(all,k); } solverSphericalGrid_CutPole = generateSpherical2DGrid_CutPole(); } /********************************************************************//** * \author Slava Merkin (vgm at bu.edu) * \since 11-2007 * * \param[in] inNLAT: Size of inLATITUDE array. * * \param[in] inLATITUDE: 1d array of latitudinal positions. * - Polar cap is at: inLatitude[0] * - Low latitude boundary is at: inLatitude[inNLAT-1] * * \param[in] LON_GRID: Longitudinal Grid Step. Requirements: * - LON_GRID <= 5 degrees * - 360 % LON_GRID ==0 and the result should be an even number. * * CREATE DOCUMENTATION LATER * ************************************************************************/ Solver::Solver(const int inNLAT, const real *inLATITUDE, const real LON_GRID) { nLat = inNLAT; nLon = 360./LON_GRID+1; positions.redim(nLat*nLon,2); x.redim(nLat,nLon); y.redim(nLat,nLon); latitude.redim(nLat,nLon); longitude.redim(nLat,nLon); realArray lon(nLon), lat(nLat); for (int i = 0; i < nLat; i++) {lat(i) = inLATITUDE[i];} lon.seqAdd(0.,LON_GRID); lon = lon/180.*PI; lat = lat/180.*PI; // ===== The part below is the same for all constructors ========== Range NLAT(0,nLat-1); for (int k = 0; k < nLon; k++) { latitude(all,k) = lat; longitude(all,k) = lon(k); x(all,k)=sin(lat)*cos(lon(k)); y(all,k)=sin(lat)*sin(lon(k)); positions(NLAT+k*nLat,0) = x(all,k); positions(NLAT+k*nLat,1) = y(all,k); } solverSphericalGrid_CutPole = generateSpherical2DGrid_CutPole(); } /********************************************************************//** * \author Slava Merkin (vgm at bu.edu) * \since 06-2007 * * \return The Solver grid projection onto XY plane * * The function creates projection of the Solver grid onto XY * plane. In this plane the interpolation between the Solver grid * functions and the MHD grid functions occurs. ************************************************************************/ CompositeGrid Solver::generateCartesian2DGrid() const { realArray ionArray(nLat,nLon,2); ionArray(all,all,0)=x; ionArray(all,all,1)=y; DataPointMapping & sphere = *new DataPointMapping; sphere.incrementReferenceCount(); sphere.setDataPoints(ionArray,2,2); sphere.setIsPeriodic(1,Mapping::functionPeriodic); CompositeGrid sgrid(2,1); sgrid[0].reference(sphere); sgrid.update(); sphere.decrementReferenceCount(); return sgrid; }; /********************************************************************//** * \author Slava Merkin (vgm at bu.edu) * \since 06-2007 * * \return The Solver grid in spherical angular coordinates * * The function creates the rectangular Solver grid in spherical * angular coordinates. This is the grid on which the potential * solution is obtained. ************************************************************************/ MappedGrid Solver::generateSpherical2DGrid() const { // SquareMapping & square = *new SquareMapping; // square.incrementReferenceCount(); // square.setVertices(0., lowLatBound/180.*PI, 0., 2*PI); // square.setGridDimensions(0,nLat); // square.setGridDimensions(1,nLon); // square.setIsPeriodic(1,Mapping::derivativePeriodic); // MappedGrid sgrid = square; // sgrid.update(); // square.decrementReferenceCount(); realArray ionArray(nLat,nLon,2); ionArray(all,all,0)=latitude; ionArray(all,all,1)=longitude; DataPointMapping & sphere = *new DataPointMapping; sphere.incrementReferenceCount(); sphere.setDataPoints(ionArray,2,2); sphere.setIsPeriodic(1,Mapping::derivativePeriodic); MappedGrid sgrid(2); sgrid.reference(sphere); sgrid.update(); sphere.decrementReferenceCount(); return sgrid; }; /********************************************************************//** * \author Slava Merkin (vgm at bu.edu) * \since 09-02-2007 * * \return The Solver grid in spherical angular coordinates but * without pole * * The function creates the rectangular Solver grid in spherical * angular coordinates but cuts out the pole. This is the grid on * which the potential solution is obtained. ************************************************************************/ MappedGrid Solver::generateSpherical2DGrid_CutPole() const { // SquareMapping & square = *new SquareMapping; // square.incrementReferenceCount(); // square.setVertices(latGrid/180.*PI, lowLatBound/180.*PI, 0., 2*PI); // square.setGridDimensions(0,nLat-1); // square.setGridDimensions(1,nLon); // square.setIsPeriodic(1,Mapping::derivativePeriodic); // MappedGrid sgrid = square; // sgrid.update(); // square.decrementReferenceCount(); // return sgrid; Range NLATm1(1,nLat-1); realArray ionArray(nLat-1,nLon,2); ionArray(all,all,0)=latitude(NLATm1,all); ionArray(all,all,1)=longitude(NLATm1,all); DataPointMapping & sphere = *new DataPointMapping; sphere.incrementReferenceCount(); sphere.setDataPoints(ionArray,2,2); sphere.setIsPeriodic(1,Mapping::derivativePeriodic); MappedGrid sgrid(2); sgrid.reference(sphere); #ifdef LOW_LAT_NEUMANN // This is needed for the Neumann BC at the low lat boundary for an // extra equation which fixes the potential solution (removes // singularity from the matrix). We set the number of ghostlines at // the lowlat boundary wall to 2 and use just one point from that // line for the extra equation. Other points remain unused sgrid.setNumberOfGhostPoints(1,0,2); #endif sgrid.update(); sphere.decrementReferenceCount(); return sgrid; }; /********************************************************************//** * \author Slava Merkin (vgm at bu.edu) * \since 06-2007 * * \param[in] fromFunc Grid function from which to interolate * * \param[in] toFunc Grid function to which to interpolate * * The two grid functions passed as parameters are assumed to be * defined on an XY grid. The toFunc should be defined on a grid * returned by generateSpherical2DGrid(), because the array #positions * created at the time of the grid generation is used for * interpolation. The fromFunc can be whatever (as long as its grid is * in XY plane), but we use only MHD grid functions so far, * i.e. MHD_Grids::getFirstShellIonCartesianGrid() or * MHD_Grids::getSecondShellIonCartesianGrid(). ************************************************************************/ void Solver::mapToSolverGrid(realCompositeGridFunction & fromFunc, realCompositeGridFunction & toFunc) { // Note, we can only interpolate contiguous components in the fromFunc this way Range R(SOLVER_FAC_NORTH,SOLVER_DENSITY_SOUTH); // Interpolate only these functions realCompositeGridFunction v; v.link(fromFunc,R); realArray u(nLat*nLon,6); u = 0.; int r = interpolatePoints(positions,v,u); Range NLAT(0,nLat-1); for (int k = 0; k < nLon; k++) { // Note, the components inside toFunc (LHS) are numbered as // functions.h describes. However, we subtract 2 from each // component on the RHS since "u" is a realArray with 6 // components. toFunc[0](NLAT,k,all,SOLVER_FAC_NORTH) = u(NLAT+k*nLat,MHD_FAC_NORTH - 2); toFunc[0](NLAT,k,all,SOLVER_FAC_SOUTH) = u(NLAT+k*nLat,MHD_FAC_SOUTH - 2); toFunc[0](NLAT,k,all,SOLVER_SOUND_SPEED_NORTH) = u(NLAT+k*nLat,MHD_SOUND_SPEED_NORTH - 2); toFunc[0](NLAT,k,all,SOLVER_SOUND_SPEED_SOUTH) = u(NLAT+k*nLat,MHD_SOUND_SPEED_SOUTH - 2); toFunc[0](NLAT,k,all,SOLVER_DENSITY_NORTH) = u(NLAT+k*nLat,MHD_DENSITY_NORTH - 2); toFunc[0](NLAT,k,all,SOLVER_DENSITY_SOUTH) = u(NLAT+k*nLat,MHD_DENSITY_SOUTH - 2); } }; /********************************************************************//** * \author Slava Merkin (vgm at bu.edu) * \since 06-2007 * * \param[in] fromFunc Grid function from which to interpolate * * \param[in] toFunc Grid function to which to interpolate * * This is the counterpart to mapToSolverGrid. Here we interpolate * from the solver grid to another (MHD) grid. This function is a * little bit more complicated than mapToSolverGrid, because we need * to extract the dimensions and positions from the grid to which we * are interpolating. ************************************************************************/ void Solver::mapFromSolverGrid(realCompositeGridFunction & fromFunc, realCompositeGridFunction & toFunc) { MappedGrid *grid = toFunc[0].getMappedGrid(); Index I1,I2,I3; getIndex((*grid).gridIndexRange(),I1,I2,I3); const realArray & x = (*grid).vertex()(I1,I2,I3,0); const realArray & y = (*grid).vertex()(I1,I2,I3,1); const int nx = x.getBound(0) - x.getBase(0) + 1; const int ny = x.getBound(1) - x.getBase(1) + 1; // OK, got the necessary info Range R(SOLVER_POTENTIAL_NORTH,SOLVER_POTENTIAL_SOUTH); // Interpolate only these functions realCompositeGridFunction v; v.link(fromFunc,R); realArray u(nx*ny,2), pos(nx*ny,2); u = 0.; // Now, create the positions array for (int k = 0; k < ny; k++) { pos(I1+k*nx,0) = x(all,k); pos(I1+k*nx,1) = y(all,k); } int r = interpolatePoints(pos,v,u); for (int k = 0; k < ny; k++) { toFunc[0](I1,k,all,MHD_POTENTIAL_NORTH) = u(I1+k*nx,SOLVER_POTENTIAL_NORTH); toFunc[0](I1,k,all,MHD_POTENTIAL_SOUTH) = u(I1+k*nx,SOLVER_POTENTIAL_SOUTH); } }; /********************************************************************//** * \author Slava Merkin (vgm at bu.edu) * \since 06-2007 * * \param[in] solverFunctions Solver grid functions defined in * functions.h * * \param[in] HEMISPHERE Should be NORTH or SOUTH (see const.h) * * The main function of the class. It obtains the solution of the for * the ionospheric potential by using Overture * OGES. We solve the equation in spherical coordinates. The * equation is [Goodman, 1995]: * * \f[ * \nabla_\perp \left[ \begin{array}{ccc} * \Sigma_P & -\Sigma_H\cos\delta \\ * \Sigma_H\cos^{-1}\delta & \Sigma_P \end{array}\right] * \nabla_\perp\Psi = j_\parallel\cos\delta * \f] * * Here \f$\Psi\f$ is the electrostatic potential, \f$j_\parallel\f$ * is the field-aligned current, \f$\Sigma_P,\,\Sigma_H\f$- Pedersen * and Hall conductances, and \f$\delta\f$ is the magnetic field dip * angle: * * \f[ * \cos\delta = -2\frac{\cos\theta}{\sqrt{1+3\cos^2\theta}} * \f] * * for the northern hemisphere, where \f$\theta\f$ is the polar angle * (colatitude). \f$\nabla_\perp\f$ is the transverse component of * divergence and gradient in spherical coordinates in the first and * second case, respectively. * * After expansion the equation above becomes * * \f[ * \frac{1}{\sin\theta}\frac{\partial}{\partial\theta}\left[ * \sin\theta\Sigma_P\frac{\partial\Psi}{\partial\theta}- * \cos\delta\Sigma_H\frac{\partial\Psi}{\partial\phi}\right]+ * \frac{\partial}{\partial\phi}\left[ * \frac{\Sigma_P}{\sin^2\theta}\frac{\partial\Psi}{\partial\phi}+ * \frac{\Sigma_H}{\sin\theta \cos\delta}\frac{\partial\Psi}{\partial\theta}\right] = * j_\parallel R^2 \cos\delta, * \f] * * where \f$R\f$ is the radius of the ionosphere. The equation is, of * course, singular at the pole. Thus the grid that is actually used * for the solution is one with the pole cut out (see * generateSpherical2DGrid_CutPole()). The grid is a square in * (\f$\theta,\phi\f$) coordinates. * * The boundary condition at the low-latitude boundary is Dirichlet: * \f$ \Psi|_{\mathrm{low \, latitude}} = 0, \f$ but it can be changed * to just about anything, e.g. \f$ * \frac{\partial\Psi}{\partial\theta} = 0\f$. * * The boundary condition at the pole is a little more sophisticated * and it has to do with the singularity of the equation. In order to * obtain a rigorous solution one needs to expand the equation above * near the pole. The resulting equation is as follows: * * \f[ * 4\sigma_P^{(0)}\Psi^{(2)} + 2\sigma_H^{(0)} \frac{\partial * \Psi^{(2)}}{\partial \phi} + \frac{\partial}{\partial \phi} * \left[\sigma_P^{(0)}\frac{\partial \Psi^{(2)}}{\partial * \phi}\right] + 2\frac{\partial}{\partial * \phi}\left[\sigma_H^{(0)}\Psi^{(2)}\right] = j_\parallel^{(0)}R^2, * \f] * * where \f$\sigma_P^{(0)}=\Sigma_P(0,\phi)$; * $\sigma_H^{(0)}=\Sigma_H(0,\phi)$; $\Psi^{(2)}=\frac{\partial^2 * \Psi}{\partial \theta^2} |_{\theta=0}\f$, and \f$\frac{\partial * \Psi}{\partial \theta}|_{\theta=0}\f$ is required for \f$\Psi\f$ to * be finite at the pole. * * In order to obtain an exact analytic solution, one would have to * integrate the equation above and match its solution with the * original equation at the pole. We simplify this task by requiring * \f$ \partial/\partial\phi|_{\theta=0}=0\f$ and then one obtains * * \f[ * 4\sigma_P^{(0)}\Psi^{(2)}=j_\parallel^{(0)}R^2, * \f] * * and therefore * * \f[ * \Psi^{(2)}=\frac{j_\parallel^{(0)}R^2}{4\sigma_P^{(0)}}. * \f] * * Thus, for a small \f$\delta\theta\f$ we have azimuthally symmetric * * \f[ * \Psi(\delta\theta) = \Psi(0,0) + \Psi^{(2)}\delta\theta^2 = * \Psi(0,0) + \frac{j_\parallel^{(0)}R^2}{4\sigma_P^{(0)}}\delta\theta^2. * \f] * * This gives us a Neumann boundary condition at a small \f$\delta\theta\f$ * boundary (circle around the pole) for the original equation: * * \f[ \frac{\partial\Psi}{\partial\theta} = * \frac{j_\parallel^{(0)}R^2}{2\sigma_P^{(0)}}\delta\theta. \f] * * And then, once the eqution is solved upto a small colatitude * \f$\delta\theta\f$, one can infer the value of the potential at the * pole: * * \f[ * \Psi(0,0) = \Psi(\delta\theta) - \frac{j_\parallel^{(0)}R^2}{2\sigma_P^{(0)}} * \delta\theta^2 * \f] * * \note As mentioned in MHD_Interface documentation, in this code, * y-axis in the southern hemisphere is inverted with respect to the * conventional LFM ionosphere, so that in the southern hemisphere * {x,y,z}-triad has the right orientation when z-axis is pointing * from the Earth's center. This way, we can apply the same solver * code in both hemispheres by just changing the sign of the dip angle. * * \b Refereces: * * M. L. Goodman, A three-dimensional, iterative mapping procedure for * the implementation of an ionosphere-magnetosphere anisotropic Ohm's * law boundary condition in global magnetohydrodynamic simulations, * Ann. Geophys., v. 13, 8, pp. 843--853, 1995. * ************************************************************************/ void Solver::solve(realCompositeGridFunction & solverFunctions, const bool HEMISPHERE) { int POTENTIAL, JPARA, SIGMAP, SIGMAH; if (HEMISPHERE==NORTH) { POTENTIAL = SOLVER_POTENTIAL_NORTH; JPARA = SOLVER_FAC_NORTH; SIGMAP = SOLVER_SIGMAP_NORTH; SIGMAH = SOLVER_SIGMAH_NORTH; } else { POTENTIAL = SOLVER_POTENTIAL_SOUTH; JPARA = SOLVER_FAC_SOUTH; SIGMAP = SOLVER_SIGMAP_SOUTH; SIGMAH = SOLVER_SIGMAH_SOUTH; } real jpara0 = solverFunctions[0](0,0,0,JPARA); real sigmap0 = solverFunctions[0](0,0,0,SIGMAP); realMappedGridFunction psi(solverSphericalGrid_CutPole), jpara(solverSphericalGrid_CutPole), sigmaP(solverSphericalGrid_CutPole), sigmaH(solverSphericalGrid_CutPole), theta(solverSphericalGrid_CutPole), cosDelta(solverSphericalGrid_CutPole), cosDelta_sinTheta(solverSphericalGrid_CutPole), RHS(solverSphericalGrid_CutPole), checkSolution(solverSphericalGrid_CutPole), // and some helper functions to hold coeffs scalar1(solverSphericalGrid_CutPole), scalar2(solverSphericalGrid_CutPole), scalar3(solverSphericalGrid_CutPole), scalar4(solverSphericalGrid_CutPole); realMappedGridFunction & vertex = solverSphericalGrid_CutPole.vertex(); theta = vertex(all,all,all,0); // This copies the data part of grid functions // As far as I understand, I cannot do an alias here as I did before (see ion_solver.C) // Because now psi and psiComp, etc. are constructed on different grids Range LAT(0,nLat); psi = 0.; // psi(LAT-1,all) = solverFunctions[0](LAT,all,all,POTENTIAL); jpara(LAT-1,all) = solverFunctions[0](LAT,all,all,JPARA); sigmaP(LAT-1,all) = solverFunctions[0](LAT,all,all,SIGMAP); sigmaH(LAT-1,all) = solverFunctions[0](LAT,all,all,SIGMAH); #ifdef TEST_SOLUTION Range J1(15,20), J2(30,35), J(25,35); real j0=1.e-6; jpara = 0.; jpara(J1,all) = j0; jpara(J2,all) = -j0; #endif // Define the dip. angle (Eq. (21), (22) Goodman (1995)) if ( HEMISPHERE == NORTH) { cosDelta = -2.*cos(theta)/sqrt( 1. + 3.*cos(theta)*cos(theta) );} else { cosDelta = 2.*cos(theta)/sqrt( 1. + 3.*cos(theta)*cos(theta) );} // We will also need this: cosDelta_sinTheta = cosDelta*sin(theta); //------------- grid functions to hold the coeffs --------------- int stencilSize=int(pow(3,solverSphericalGrid_CutPole.numberOfDimensions())); realMappedGridFunction coeff(solverSphericalGrid_CutPole,stencilSize,all,all,all), coeff1(solverSphericalGrid_CutPole,stencilSize,all,all,all), coeff2(solverSphericalGrid_CutPole,stencilSize,all,all,all), coeff3(solverSphericalGrid_CutPole,stencilSize,all,all,all), coeff4(solverSphericalGrid_CutPole,stencilSize,all,all,all); coeff.setIsACoefficientMatrix(TRUE,stencilSize); coeff1.setIsACoefficientMatrix(TRUE,stencilSize); coeff2.setIsACoefficientMatrix(TRUE,stencilSize); coeff3.setIsACoefficientMatrix(TRUE,stencilSize); coeff4.setIsACoefficientMatrix(TRUE,stencilSize); //---------------------------------------------------------------- //------------ mapped grid operators ---------------------------- MappedGridOperators op(solverSphericalGrid_CutPole); op.setStencilSize(stencilSize); coeff.setOperators(op); coeff1.setOperators(op); coeff2.setOperators(op); coeff3.setOperators(op); coeff4.setOperators(op); scalar1.setOperators(op); scalar2.setOperators(op); scalar3.setOperators(op); scalar4.setOperators(op); cosDelta_sinTheta.setOperators(op); RHS.setOperators(op); //---------------------------------------------------------------- //------------ equation ------------------------------------------ const int pole = 1, lowLatitudeBoundary = 2; /* side axis */ solverSphericalGrid_CutPole.boundaryCondition()( 0, 0) = pole; solverSphericalGrid_CutPole.boundaryCondition()( 1, 0) = lowLatitudeBoundary; // VGM 06-13-08. It is important to extrapolate the ghostline on the // low-latitude boundary (applying BCs to the coefficient matrix is // not enough). I've found this out after I realized that the // solution I obtained for Neuman BC was not forced exactly. This // fixes the problem sigmaP.setOperators(op); sigmaH.setOperators(op); sigmaP.applyBoundaryCondition(0,BCTypes::extrapolate,lowLatitudeBoundary); sigmaP.finishBoundaryConditions(); sigmaH.applyBoundaryCondition(0,BCTypes::extrapolate,lowLatitudeBoundary); sigmaH.finishBoundaryConditions(); /* Periodic boundary is set up by default from the mapping - no need to set it here manually */ // scalar1 = sigmaP*sin(theta); // scalar2 = -sigmaH*cosDelta; // scalar3 = sigmaH/cosDelta/sin(theta); // scalar4 = sigmaP/sin(theta)/sin(theta); scalar1 = sigmaP*sin(theta)/cosDelta/cosDelta; scalar2 = -sigmaH/cosDelta; scalar3 = sigmaH*sin(theta); scalar4 = sigmaP*cosDelta; scalar1.finishBoundaryConditions(); scalar2.finishBoundaryConditions(); scalar3.finishBoundaryConditions(); scalar4.finishBoundaryConditions(); coeff1 = op.derivativeScalarDerivativeCoefficients(scalar1,0,0); coeff2 = op.derivativeScalarDerivativeCoefficients(scalar2,0,1); coeff3 = op.derivativeScalarDerivativeCoefficients(scalar3,1,0); coeff4 = op.derivativeScalarDerivativeCoefficients(scalar4,1,1); coeff = multiply(cosDelta_sinTheta,coeff1) + multiply(cosDelta_sinTheta,coeff2) + coeff3 + coeff4; RHS = jpara*pow(cosDelta_sinTheta*rIonosphere,2); // coeff = multiply(1./sin(theta),coeff1) + // multiply(1./sin(theta),coeff2) + // coeff3 + // coeff4; // RHS = jpara*rIonosphere*rIonosphere*cosDelta; #ifndef LOW_LAT_NEUMANN coeff.applyBoundaryConditionCoefficients(0,0,BCTypes::dirichlet,lowLatitudeBoundary); coeff.applyBoundaryConditionCoefficients(0,0,BCTypes::extrapolate,lowLatitudeBoundary); #else // This is for Neumann at low lat boundary coeff.applyBoundaryConditionCoefficients(0,0,BCTypes::neumann,lowLatitudeBoundary); #endif coeff.applyBoundaryConditionCoefficients(0,0,BCTypes::neumann,pole); coeff.finishBoundaryConditions(); // Periodic RHS for ( int i = 0; i < solverSphericalGrid_CutPole.numberOfGhostPoints(0,1) + 1; i++) { RHS(all,-i) = 0.; } for ( int i = 1; i < solverSphericalGrid_CutPole.numberOfGhostPoints(1,1) + 1; i++) { RHS(all,nLon-1+i) = 0.; } Index Ib1, Ib2, Ib3; Index Ig1, Ig2, Ig3; // Low lat. RHS #ifndef LOW_LAT_NEUMANN getBoundaryIndex(solverSphericalGrid_CutPole.gridIndexRange(),End,axis1,Ib1,Ib2,Ib3); RHS(Ib1,Ib2,Ib3) = 0.; /* or any other boundary value for the potential */ #else getGhostIndex(solverSphericalGrid_CutPole.gridIndexRange(),End,axis1,Ig1,Ig2,Ig3); RHS(Ig1,Ig2,Ig3) = 0.; // This is for Neumann #endif // Pole RHS getGhostIndex(solverSphericalGrid_CutPole.gridIndexRange(),Start,axis1,Ig1,Ig2,Ig3); if (HEMISPHERE==NORTH) { RHS(Ig1,Ig2,Ig3) = -jpara0*rIonosphere*rIonosphere/(2.*sigmap0)*theta(0,0); } else { RHS(Ig1,Ig2,Ig3) = jpara0*rIonosphere*rIonosphere/(2.*sigmap0)*theta(0,0); } //-------------- Solver ------------------------------------------ Oges solver(solverSphericalGrid_CutPole); solver.set(OgesParameters::THEsolverType,OgesParameters::SLAP); // solver.set(OgesParameters::THEsolverMethod,OgesParameters::generalizedMinimalResidual); // solver.set(OgesParameters::THEpreconditioner,OgesParameters::incompleteLUPreconditioner); // solver.set(OgesParameters::THErelativeTolerance,1.e-11); // solver.set(OgesParameters::THEmaximumNumberOfIterations,10000); // solver.set(OgesParameters::THEnumberOfIncompleteLULevels,iluLevels); solver.setCoefficientArray(coeff); // solver.set(OgesParameters::THEabsoluteTolerance,1.e-6); #ifdef LOW_LAT_NEUMANN // We fix the singular matrix here solver.set(OgesParameters::THEcompatibilityConstraint,TRUE); // Tell the solver to refactor the matrix since the coefficients have changed solver.setRefactor(TRUE); // we need to reorder too because the matrix changes a lot for the singular case solver.setReorder(TRUE); solver.initialize(); int ne,i1e,i2e,i3e,gride; solver.equationToIndex( solver.extraEquationNumber(0),ne,i1e,i2e,i3e,gride); RHS(i1e,i2e,i3e)=0.; #endif solver.solve(psi,RHS); //---------------------------------------------------------------- solverFunctions[0](LAT,all,all,POTENTIAL) = psi(LAT-1,all); if (HEMISPHERE==NORTH) { solverFunctions[0](0,all,all,POTENTIAL) = sum(psi(0,all))/nLon + jpara0*rIonosphere*rIonosphere/(4.*sigmap0)*theta(0,0)*theta(0,0); } else { solverFunctions[0](0,all,all,POTENTIAL) = sum(psi(0,all))/nLon - jpara0*rIonosphere*rIonosphere/(4.*sigmap0)*theta(0,0)*theta(0,0); } /* psi.setOperators(op); psi.applyBoundaryCondition(0,BCTypes::dirichlet,lowLatitudeBoundary,0.); psi.applyBoundaryCondition(0,BCTypes::extrapolate,lowLatitudeBoundary,0.,0.,extrapParams); psi.applyBoundaryCondition(0,BCTypes::neumann,pole,0.); psi.finishBoundaryConditions(); checkSolution = sin(theta)*sin(theta)*psi.xx()+sin(theta)*cos(theta)*psi.x()+ psi.yy() - jpara*rIonosphere*rIonosphere/sigmaP*cosDelta*sin(theta)*sin(theta); */ #ifdef TEST_SOLUTION ofstream file("check_solution.dat",ios::out); for (int j = 0; j < nLat-1; j++) { file << theta(j,0) << "\t" << jpara(j,0) << "\t" << psi(j,0) << "\n"; } file.close(); #endif } /********************************************************************//** * \author Slava Merkin (vgm at bu.edu) * \since 06-2007 * * \return Creates a full-globe grid in SM spherical angular * coordinates * * In order to interpolate between the solver grid and the ITM grid we * need to cover the entire globe. This way, the interpolation can be * done in 2D plane of spherical angular coordinates. The angular * coordinates are defined here as in ITM_Grids rather than the * regular spherical coordinates used to solve(), i.e. latitude goes * from \f$-\pi/2\f$ (south pole) to \f$\pi/2\f$ (north pole) * longitude goes from 0 (noon) to \f$\pi=-\pi\f$ (midnight) ************************************************************************/ CompositeGrid Solver::generateSMGlobeGrid() { real deltam = latitude(nLat-1,0) - latitude(nLat-2,0); real thetam = latitude(nLat-1,0); int N = floor( (PI-2*thetam)/deltam ); real delta = (PI-2*thetam)/N; nLatGlobe = 2*nLat+ N - 1; // OK, now that we have the dimensions of the grid - let's fill in // the vertices realArray smGlobeLatitude(nLatGlobe), lattmp(N); lattmp.seqAdd(thetam+delta,delta); Range NLAT(0,nLat-1), NRANGE(nLat,nLat+N-1); smGlobeLatitude(NLAT) = latitude(NLAT,0); smGlobeLatitude(NRANGE) = lattmp; for (int k = 0; k < nLat; k++) { smGlobeLatitude(nLatGlobe-nLat+k) = PI - latitude(nLat-1-k,0); } realArray ionArray(nLatGlobe,nLon,2); for (int k = 0; k < nLon; k++) { ionArray(all,k,0) = smGlobeLatitude-PI/2.; ionArray(all,k,1) = longitude(0,k)-PI; } DataPointMapping & sphere = *new DataPointMapping; sphere.incrementReferenceCount(); sphere.setDataPoints(ionArray,2,2); sphere.setIsPeriodic(1,Mapping::derivativePeriodic); CompositeGrid sgrid(2,1); sgrid[0].reference(sphere); sgrid.update(); sphere.decrementReferenceCount(); return sgrid; }; /********************************************************************//** * \author Slava Merkin (vgm at bu.edu) * \since 06-2007 * * \param[in] fNorth A grid function on Solver grid for northern * hemisphere * * \param[in] fSouth A grid function on Solver grid for southern * hemisphere * * \param[out] fGlobe A grid function on a globe * grid=generateSMGlobeGrid() * * This function takes as inputs grid functions for north and south * defined on solver grid and fills in values over entire globe * (values at midlatitudes, in between the two polar caps are set to * 0). * * \note The indexing is not trivial because the globe grid coming out * of generateSMGlobeGrid() defines angular coordinates differently * than the Solver grid. ************************************************************************/ void Solver::mapToSMGlobe(realCompositeGridFunction & fSolver, realCompositeGridFunction & fGlobe) const { // Initialize to zero fGlobe = 0.; Range NLAT(0,nLat-1), NLON2(0,(nLon-1)/2), NLON22 ((nLon-1)/2,nLon-1); // DO THE NORTH for (int i=0; i