#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 vmg.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 vmg.c. Referenced by mgsolve(), and mgvcycle(). |
|
||||||||||||
|
Definition at line 42 of file vmg.c. References MIN. Referenced by initialize().
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 vmg.c. References mymalloc(), real, and real2D. Referenced by initialize().
00070 {
00071 int j;
00072 int err = 0;
00073
00074 err |= ((irhs = (real2D *) malloc((ng + 1) * sizeof(real2D))) == NULL);
00075 err |= ((ip = (real2D *) malloc((ng + 1) * sizeof(real2D))) == NULL);
00076 err |= ((iap = (real2D *) malloc((ng + 1) * sizeof(real2D))) == NULL);
00077 err |= ((ic0 = (real2D *) malloc((ng + 1) * sizeof(real2D))) == NULL);
00078 err |= ((ic1 = (real2D *) malloc((ng + 1) * sizeof(real2D))) == NULL);
00079 err |= ((ic2 = (real2D *) malloc((ng + 1) * sizeof(real2D))) == NULL);
00080 err |= ((ic3 = (real2D *) malloc((ng + 1) * sizeof(real2D))) == NULL);
00081 err |= ((ic4 = (real2D *) malloc((ng + 1) * sizeof(real2D))) == NULL);
00082 err |= ((ierr = (real2D *) malloc((ng + 1) * sizeof(real2D))) == NULL);
00083
00084 if ( err ) return MEMERROR;
00085
00086 err |= ((sxp = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL);
00087 err |= ((syp = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL);
00088
00089 err |= ((syu = (real2D)mymalloc(sizeof(real), n1 + 1, n2)) == NULL);
00090 err |= ((au = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL);
00091 err |= ((cu = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL);
00092
00093 err |= ((sxv = (real2D)mymalloc(sizeof(real), n1, n2 + 1)) == NULL);
00094 err |= ((av = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL);
00095 err |= ((cv = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL);
00096
00097 err |= ((rz1dr = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL);
00098 err |= ((rz2dr = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL);
00099 err |= ((rr1dz = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL);
00100 err |= ((rr2dz = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL);
00101
00102 err |= ((w1 = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL);
00103 err |= ((w2 = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL);
00104 err |= ((dummy = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL);
00105
00106 err |= ((ierr[ng] = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL);
00107 err |= ((irhs[ng] = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL);
00108
00109 n1 = n1 / 2 + 1;
00110 n2 = n2 / 2 + 1;
00111
00112 for (j = ng - 1; j >= 1; j--) {
00113 err |= ((irhs[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 err |= ((ierr[j] = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL);
00122
00123 if ( err ) return MEMERROR;
00124
00125 n1 = n1 / 2 + 1;
00126 n2 = n2 / 2 + 1;
00127 }
00128
00129 return 0;
00130 } /* end initpressure() */
|
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|
Definition at line 134 of file vmg.c. References add(), interface::bc, bc_scalar(), COARSEST, coli(), copy(), cut_interface(), fill0(), 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(). Referenced by pressure().
00138 {
00139 int ngrid, nn1, nn2, i, j;
00140 interface intemp;
00141 interface_cut cutp, cutu, cutv;
00142 real scale = 1.0;
00143 int ncycle = 0;
00144
00145 resid(irhs[ng], p, u, ap, c0, c1, c2, c3, c4, n1, n2);
00146 bc_scalar(irhs[ng], n1, n2, NUL);
00147 if (rdivmax(irhs[ng], n1, n2, &i, &j) <= maxdiv)
00148 return 0;
00149
00150 ip[ng] = p;
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 copy(w1, p, n1, n2);
00212 while (rdivmax(irhs[ng], n1, n2, &i, &j) > maxdiv && ncycle < NCYCLEMAX) {
00213 fill0(ip[ng], n1, n2);
00214 mgvcycle(ng, n1, n2, 1, 1);
00215 printf("ncycle: %d rdivmax: %g\n",
00216 ncycle + 1, rdivmax(irhs[ng], n1, n2, &i, &j));
00217 add(w1, p, n1, n2);
00218 ncycle++;
00219 }
00220 copy(p, w1, n1, n2);
00221
00222 return ncycle;
00223 } /* end mginit() */
|
|
||||||||||||||||||||||||
|
Definition at line 227 of file vmg.c. References add(), bc_scalar(), COARSEST, fill0(), NUL, NULGRAD, NULGRAD2, relax(), and resid(). Referenced by mgsolve(), and mgvcycle().
00228 {
00229 int count;
00230 int ncx = nx/2 + 1, ncy = ny/2 + 1;
00231 char s[256];
00232
00233 /* solve exactly on the coarsest grid */
00234 if (n == COARSEST) {
00235 slvsml(ip[COARSEST], irhs[COARSEST], iap[COARSEST],
00236 ic0[COARSEST], ic1[COARSEST], ic2[COARSEST], ic3[COARSEST],
00237 ic4[COARSEST],
00238 nx, ny);
00239 return;
00240 }
00241
00242 /* do prerelaxation */
00243 for (count = 0; count < npre; count++) {
00244 relax(ip[n], irhs[n], iap[n],
00245 ic0[n], ic1[n], ic2[n], ic3[n], ic4[n], nx, ny);
00246 bc_scalar(ip[n], nx, ny, NULGRAD2);
00247 }
00248
00249 /* compute defect */
00250 resid(irhs[n], ip[n], irhs[n], iap[n],
00251 ic0[n], ic1[n], ic2[n], ic3[n], ic4[n],
00252 nx, ny);
00253 bc_scalar(irhs[n], nx, ny, NUL);
00254
00255 /* compute r.h.s. on the coarser grid using defect */
00256 rstrct(irhs[n-1], irhs[n], ncx, ncy, - 4.0);
00257
00258 /* use zero as initial guess on the coarser grid */
00259 fill0(ip[n-1], ncx, ncy);
00260
00261 /* solve on the coarser grid */
00262 mgvcycle(n - 1, ncx, ncy, npre + 1, npost + 1);
00263
00264 /* interpolate the error from coarse to fine grid */
00265 interp(ierr[n], ip[n-1], nx, ny);
00266
00267 /* new residual res(new) = res(old)-A*error */
00268 resid(irhs[n], ierr[n], irhs[n], iap[n],
00269 ic0[n], ic1[n], ic2[n], ic3[n], ic4[n],
00270 nx, ny);
00271 bc_scalar(irhs[n], nx, ny, NUL);
00272
00273 /* correct e(ngrid) with error = interp(e(ngrid-1)) */
00274 add(ip[n], ierr[n], nx, ny);
00275
00276 /* do post relaxation with initial guess = 0 */
00277 fill0(ierr[n], nx, ny);
00278 for (count = 0; count < npost; count++) {
00279 relax(ierr[n], irhs[n], iap[n],
00280 ic0[n], ic1[n], ic2[n], ic3[n], ic4[n], nx, ny);
00281 bc_scalar(ierr[n], nx, ny, NULGRAD);
00282 }
00283
00284 /* add correction to correction to solution */
00285 add(ip[n], ierr[n], nx, ny);
00286
00287 /* compute residual */
00288 resid(irhs[n], ierr[n], irhs[n], iap[n],
00289 ic0[n], ic1[n], ic2[n], ic3[n], ic4[n],
00290 nx, ny);
00291 bc_scalar(irhs[n], nx, ny, NUL);
00292 /*
00293 sprintf(s, "ip%d", n);
00294 datbarray(ip[n], nx, ny, s);
00295 sprintf(s, "ipc%d", n);
00296 sprintf(s, "iap%d", n);
00297 datbarray(iap[n], nx, ny, s);
00298 */
00299 }
|
|
||||||||||||||||||||||||||||||||||||||||||||
|
Definition at line 409 of file vmg.c. Referenced by mgfas(), mgvcycle(), and pressure().
00412 {
00413 int i, j, isw, jsw;
00414
00415 for (jsw = 1; jsw <= 2; jsw++) {
00416 isw = jsw;
00417 for (i = 2; i < nx; i++, isw = 3 - isw)
00418 for (j = isw + 1; j < ny; j += 2)
00419 if (INP(i,j))
00420 p[i][j] = (c1[i][j]*p[i+1][j] + c2[i][j]*p[i][j+1] +
00421 c3[i][j]*p[i-1][j] + c4[i][j]*p[i][j-1]
00422 - rhs[i][j])/c0[i][j];
00423 }
00424 }
|
|
||||||||||||||||||||||||||||||||||||||||||||||||
|
Definition at line 391 of file vmg.c. Referenced by mgsolve(), and mgvcycle().
00394 {
00395 int i, j;
00396
00397 for (i = 2; i < nx; i++)
00398 for (j = 2; j < ny; j++)
00399 if (INP(i,j))
00400 res[i][j] = c1[i][j]*p[i+1][j] + c2[i][j]*p[i][j+1]
00401 + c3[i][j]*p[i-1][j] + c4[i][j]*p[i][j-1] - c0[i][j]*p[i][j]
00402 - rhs[i][j];
00403 else
00404 res[i][j] = 0.0;
00405 }
|
1.2.18