#include "utilf.h"#include "mymalloc.h"#include "markers.h"#include "surfaces.h"#include "extra.h"#include "make_bc.h"#include "momentum.h"#include "mg.h"#include "inout.h"Include dependency graph for mg.c:
Go to the source code of this file.
Defines | |
| #define | COARSEST 3 |
| #define | COARSERELAX 20 |
Functions | |
| int | computeng (int n1, int n2) |
| int | initpressure (int n1, int n2, int ng) |
| int | mgsolve (real2D u, real2D p, real2D ap, real2D c0, real2D c1, real2D c2, real2D c3, real2D c4, interface in, real lredis, real maxdiv, int n1, int n2, int ng) |
| void | mgvcycle (int n, int nx, int ny, int npre, int npost) |
| void | resid (real2D res, real2D p, real2D rhs, real2D ap, real2D c0, real2D c1, real2D c2, real2D c3, real2D c4, int nx, int ny) |
| void | relax (real2D p, real2D rhs, real2D ap, real2D c0, real2D c1, real2D c2, real2D c3, real2D c4, int nx, int ny) |
|
|
|
|
|
Definition at line 19 of file mg.c. Referenced by mgsolve(), and mgvcycle(). |
|
||||||||||||
|
Definition at line 42 of file mg.c. References MIN.
00043 {
00044 int nmin, n, nn, j;
00045
00046 nmin = MIN(n1, n2);
00047
00048 j = nmin - 2; n = 1;
00049 while (j > 2)
00050 { j /= 2; n++; }
00051
00052 nn = 1;
00053 for (j = 1; j <= n; j++)
00054 nn *= 2;
00055 if((nn + 2) != nmin)
00056 return NGERROR;
00057
00058 if ((n1 - 2) % (nmin - 2) != 0)
00059 return NXERROR;
00060 if ((n2 - 2) % (nmin - 2) != 0)
00061 return NYERROR;
00062
00063 Ngrid = n;
00064 return n;
00065 }
|
|
||||||||||||||||
|
Definition at line 69 of file mg.c. References mymalloc(), real, and real2D.
00070 {
00071 int j;
00072 int err = 0;
00073
00074 err |= ((irhs = (real2D *) malloc((ng + 1) * sizeof(real2D))) == NULL);
00075 err |= ((ires = (real2D *) malloc((ng + 1) * sizeof(real2D))) == NULL);
00076 err |= ((ip = (real2D *) malloc((ng + 1) * sizeof(real2D))) == NULL);
00077 err |= ((iap = (real2D *) malloc((ng + 1) * sizeof(real2D))) == NULL);
00078 err |= ((ic0 = (real2D *) malloc((ng + 1) * sizeof(real2D))) == NULL);
00079 err |= ((ic1 = (real2D *) malloc((ng + 1) * sizeof(real2D))) == NULL);
00080 err |= ((ic2 = (real2D *) malloc((ng + 1) * sizeof(real2D))) == NULL);
00081 err |= ((ic3 = (real2D *) malloc((ng + 1) * sizeof(real2D))) == NULL);
00082 err |= ((ic4 = (real2D *) malloc((ng + 1) * sizeof(real2D))) == NULL);
00083
00084 if ( err ) return MEMERROR;
00085
00086 err |= ((ires[ng] = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL);
00087
00088 err |= ((sxp = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL);
00089 err |= ((syp = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL);
00090
00091 err |= ((syu = (real2D)mymalloc(sizeof(real), n1 + 1, n2)) == NULL);
00092 err |= ((au = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL);
00093 err |= ((cu = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL);
00094
00095 err |= ((sxv = (real2D)mymalloc(sizeof(real), n1, n2 + 1)) == NULL);
00096 err |= ((av = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL);
00097 err |= ((cv = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL);
00098
00099 err |= ((rz1dr = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL);
00100 err |= ((rz2dr = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL);
00101 err |= ((rr1dz = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL);
00102 err |= ((rr2dz = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL);
00103
00104 err |= ((w1 = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL);
00105 err |= ((w2 = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL);
00106 err |= ((dummy = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL);
00107
00108 n1 = n1 / 2 + 1;
00109 n2 = n2 / 2 + 1;
00110
00111 for (j = ng - 1; j >= 1; j--) {
00112 err |= ((irhs[j] = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL);
00113 err |= ((ires[j] = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL);
00114 err |= ((ip[j] = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL);
00115 err |= ((iap[j] = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL);
00116 err |= ((ic0[j] = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL);
00117 err |= ((ic1[j] = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL);
00118 err |= ((ic2[j] = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL);
00119 err |= ((ic3[j] = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL);
00120 err |= ((ic4[j] = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL);
00121
00122 if ( err ) return MEMERROR;
00123
00124 n1 = n1 / 2 + 1;
00125 n2 = n2 / 2 + 1;
00126 }
00127
00128 return 0;
00129 } /* end initpressure() */
|
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|
Definition at line 133 of file mg.c. References interface::bc, bc_scalar(), COARSEST, coli(), cut_interface(), interface_free(), interface_init(), interface_redis_scale(), interface_rrdz(), interface_rzdr(), interface_splines(), interface_surfaces(), interface_cut::list, mgvcycle(), NCYCLEMAX, NUL, rdivmax(), real, real2D, resid(), sx_fill(), and sy_fill().
00137 {
00138 int ngrid, nn1, nn2, i, j;
00139 interface intemp;
00140 interface_cut cutp, cutu, cutv;
00141 real scale = 1.0;
00142 int ncycle = 0;
00143
00144 resid(ires[ng], p, u, ap, c0, c1, c2, c3, c4, n1, n2);
00145 bc_scalar(ires[ng], n1, n2, NUL);
00146 if (rdivmax(ires[ng], n1, n2, &i, &j) <= maxdiv)
00147 return 0;
00148
00149 ip[ng] = p;
00150 irhs[ng] = u;
00151 iap[ng] = ap;
00152 ic0[ng] = c0; ic1[ng] = c1; ic2[ng] = c2; ic3[ng] = c3; ic4[ng] = c4;
00153 interface_init(&intemp);
00154 intemp.bc = in.bc;
00155 cutp.list = cutu.list = cutv.list = NULL;
00156
00157 nn1 = n1;
00158 nn2 = n2;
00159 ngrid = ng;
00160 while (ngrid > COARSEST) {
00161 nn1 = nn1 / 2 + 1;
00162 nn2 = nn2 / 2 + 1;
00163 ngrid--;
00164 scale *= 2.0;
00165 lredis *= 2.0;
00166 /* rescale the interface */
00167 interface_redis_scale(&intemp, in, lredis, 1./scale);
00168 interface_splines(intemp, nn1, nn2);
00169 /* compute volume fraction */
00170 rstrct(iap[ngrid], iap[ngrid+1], nn1, nn2, 1.0);
00171
00172 cut_interface(intemp, &cutp, 0.0, 0.0);
00173 cut_interface(intemp, &cutu, -0.5, 0.0);
00174 cut_interface(intemp, &cutv, 0.0, -0.5);
00175
00176 interface_surfaces(intemp, cutp, sxp, syp, iap[ngrid], dummy,
00177 w1, w2, nn1, nn2);
00178 for (i = 2; i <= nn1; i++)
00179 for (j = 2; j <= nn2; j++) {
00180 au[i][j] = 0.5*(iap[ngrid][i][j] + iap[ngrid][i-1][j]);
00181 av[i][j] = 0.5*(iap[ngrid][i][j] + iap[ngrid][i][j-1]);
00182 }
00183 interface_surfaces(intemp, cutu, dummy, syu, au, cu,
00184 w1, w2, nn1, nn2);
00185 interface_surfaces(intemp, cutv, sxv, dummy, av, cv,
00186 w1, w2, nn1, nn2);
00187 sx_fill(sxp, iap[ngrid], 0.0, 0.0, intemp, nn1, nn2);
00188 sy_fill(syp, iap[ngrid], 0.0, 0.0, intemp, nn1, nn2);
00189 sy_fill(syu, au, -0.5, 0.0, intemp, nn1, nn2);
00190 for (j = 2; j <= nn2; j++)
00191 syu[nn1+1][j] = (real)j - 1.5;
00192 sx_fill(sxv, av, 0.0, -0.5, intemp, nn1, nn2);
00193 for (i = 2; i <= nn1; i++)
00194 sxv[i][nn2+1] = (real)nn2 - 1.5;
00195
00196 interface_rzdr(intemp, cutp, -0.5, rz1dr, nn1, nn2);
00197 interface_rzdr(intemp, cutp, 0.5, rz2dr, nn1, nn2);
00198 interface_rrdz(intemp, cutp, -0.5, rr1dz, nn1, nn2);
00199 interface_rrdz(intemp, cutp, 0.5, rr2dz, nn1, nn2);
00200
00201 coli(sxp, syp, syu, sxv, cu, cv, iap[ngrid], av,
00202 rz1dr, rz2dr, rr1dz, rr2dz,
00203 ic0[ngrid], ic1[ngrid], ic2[ngrid], ic3[ngrid], ic4[ngrid],
00204 nn1, nn2, scale);
00205 }
00206 interface_free(intemp);
00207 if (cutp.list) free(cutp.list);
00208 if (cutu.list) free(cutu.list);
00209 if (cutv.list) free(cutv.list);
00210
00211 while (rdivmax(ires[ng], n1, n2, &i, &j) > maxdiv && ncycle < NCYCLEMAX) {
00212 mgvcycle(ng, n1, n2, 4, 4);
00213 resid(ires[ng], p, u, ap, c0, c1, c2, c3, c4, n1, n2);
00214 bc_scalar(ires[ng], n1, n2, NUL);
00215 ncycle++;
00216 }
00217 return ncycle;
00218 } /* end mgsolve() */
|
|
||||||||||||||||||||||||
|
Definition at line 222 of file mg.c. References bc_scalar(), COARSEST, fill0(), mgvcycle(), NUL, NULGRAD2, relax(), and resid().
00223 {
00224 int count;
00225 int ncx = nx/2 + 1, ncy = ny/2 + 1;
00226
00227 /* solve exactly on the coarsest grid */
00228 if (n == COARSEST) {
00229 slvsml(ip[COARSEST], irhs[COARSEST], iap[COARSEST],
00230 ic0[COARSEST], ic1[COARSEST], ic2[COARSEST], ic3[COARSEST],
00231 ic4[COARSEST],
00232 nx, ny);
00233 return;
00234 }
00235
00236 /* do prerelaxation */
00237 for (count = 0; count < npre; count++) {
00238 relax(ip[n], irhs[n], iap[n],
00239 ic0[n], ic1[n], ic2[n], ic3[n], ic4[n], nx, ny);
00240 bc_scalar(ip[n], nx, ny, NULGRAD2);
00241 }
00242
00243 /* compute defect */
00244 resid(ires[n], ip[n], irhs[n], iap[n],
00245 ic0[n], ic1[n], ic2[n], ic3[n], ic4[n],
00246 nx, ny);
00247 bc_scalar(ires[n], nx, ny, NUL);
00248
00249 /* compute r.h.s. on the coarser grid using defect */
00250 rstrct(irhs[n-1], ires[n], ncx, ncy, - 4.0);
00251
00252 /* use zero as initial guess on the coarser grid */
00253 fill0(ip[n-1], ncx, ncy);
00254
00255 /* solve on the coarser grid */
00256 mgvcycle(n - 1, ncx, ncy, npre + 1, npost + 1);
00257
00258 /* correct on the fine grid (use ires as temporary storage) */
00259 addint(ip[n], ip[n-1], ires[n], iap[n], nx, ny);
00260
00261 /* do post relaxation */
00262 for (count = 0; count < npost; count++) {
00263 relax(ip[n], irhs[n], iap[n],
00264 ic0[n], ic1[n], ic2[n], ic3[n], ic4[n], nx, ny);
00265 bc_scalar(ip[n], nx, ny, NULGRAD2);
00266 }
00267 /*
00268 sprintf(s, "ip%d", n);
00269 datbarray(ip[n], nx, ny, s);
00270 sprintf(s, "ipc%d", n);
00271 interp(ires[n], ip[n-1], nx, ny);
00272 datbarray(ires[n], nx, ny, s);
00273 sprintf(s, "iap%d", n);
00274 datbarray(iap[n], nx, ny, s);
00275 */
00276 }
|
|
||||||||||||||||||||||||||||||||||||||||||||
|
Definition at line 386 of file mg.c.
00389 {
00390 int i, j, isw, jsw;
00391
00392 for (jsw = 1; jsw <= 2; jsw++) {
00393 isw = jsw;
00394 for (i = 2; i < nx; i++, isw = 3 - isw)
00395 for (j = isw + 1; j < ny; j += 2)
00396 if (INP(i,j))
00397 p[i][j] = (c1[i][j]*p[i+1][j] + c2[i][j]*p[i][j+1] +
00398 c3[i][j]*p[i-1][j] + c4[i][j]*p[i][j-1]
00399 - rhs[i][j])/c0[i][j];
00400 }
00401 }
|
|
||||||||||||||||||||||||||||||||||||||||||||||||
|
Definition at line 368 of file mg.c.
00371 {
00372 int i, j;
00373
00374 for (i = 2; i < nx; i++)
00375 for (j = 2; j < ny; j++)
00376 if (INP(i,j))
00377 res[i][j] = c1[i][j]*p[i+1][j] + c2[i][j]*p[i][j+1]
00378 + c3[i][j]*p[i-1][j] + c4[i][j]*p[i][j-1] - c0[i][j]*p[i][j]
00379 - rhs[i][j];
00380 else
00381 res[i][j] = 0.0;
00382 }
|
1.2.18