Why mass lumping fails in 3d?

General Discussion
1

Hello, I’m testing mass lumping in 3d for DMP but it seems to not have an effect.

heres the code :

// Small 3D heat equation example (implicit Euler)
// PARAM -wg -ns

load "msh3"

verbosity = 0;

// Mesh and time parameters
int nx = 8;
int ny = 8;
int nz = 8;
real dt = 0.001;
real t = 0.0;
real tfin = 1.0;
real eps = 1e-4;
real bx = 1.0;
real by = 0.5;
real bz = 0.5;
int useLumping = 1; // set to 1 to enable mass lumping

// Unit cube mesh
mesh3 Th = cube(nx, ny, nz);

fespace Vh(Th, P1);
Vh u, v, uold;

// Initial condition: sharp Gaussian pulse
real x0 = 0.25;
real y0 = 0.25;
real z0 = 0.25;
real sigma = 0.05;
uold = exp(-((x - x0)^2 + (y - y0)^2 + (z - z0)^2) / (2.0 * sigma^2));
real bcVal = 0.0;

// Variational forms: (u - uold)/dt + b·grad(u) - eps*Laplace(u) = 0
varf heatCons(u, v) = int3d(Th)(u * v)
                   + int3d(Th)(dt * (bx * dx(u) + by * dy(u) + bz * dz(u)) * v)
                   + int3d(Th)(dt * eps * (dx(u) * dx(v) + dy(u) * dy(v) + dz(u) * dz(v)))
                   + on(1,2,3,4,5,6, u = bcVal);
varf heatLump(u, v) = int3d(Th, qfV=qfV1lump)(u * v)
                   + int3d(Th)(dt * (bx * dx(u) + by * dy(u) + bz * dz(u)) * v)
                   + int3d(Th)(dt * eps * (dx(u) * dx(v) + dy(u) * dy(v) + dz(u) * dz(v)))
                   + on(1,2,3,4,5,6, u = bcVal);
varf massRhsCons(u,v) = int3d(Th)(uold * v)+ on(1,2,3,4,5,6, u = bcVal);
varf massRhsLump(u,v) = int3d(Th, qfV=qfV1lump)(uold * v)+ on(1,2,3,4,5,6, u = bcVal);

matrix A;
if (useLumping) {
    A = heatLump(Vh, Vh);
} else {
    A = heatCons(Vh, Vh);
}
set(A, solver = sparsesolver);

real[int] b(Vh.ndof);

real initMin = uold[].min;
real initMax = uold[].max;
real boundMin = min(initMin, bcVal);
real boundMax = max(initMax, bcVal);
real tol = 1e-12;
int violations = 0;

for (int i = 1; i * dt <= tfin; i++) {
    t = i * dt;
    if (useLumping) {
        b = massRhsLump(0, Vh);
    } else {
        b = massRhsCons(0, Vh);
    }
    u[] = A^-1 * b;
    uold = u;

    if (i % 10 == 0) {
        real minv = u[].min;
        real maxv = u[].max;
        plot(u, fill = true, value = true, cmm = "Time: " + t + " | lumping: " + useLumping, wait=true);
        if (minv < boundMin - tol || maxv > boundMax + tol) {
            violations++;
        }
        cout << "Time: " << t << " | min: " << minv << " | max: " << maxv
             << " | DMP violations: " << violations << endl;
    }
}

cout << "Initial bounds: [" << initMin << ", " << initMax << "]" << endl;
cout << "Boundary value: " << bcVal << endl;
cout << "DMP bounds: [" << boundMin << ", " << boundMax << "]" << endl;
cout << "Total DMP violations: " << violations << endl;

and heres the results :

FreeFem++.exe heat3d_small.edp -cd -wg -ns
– FreeFem++ v4.15 (Wed Dec 11 16:48:56 CET 2024 - git v4.15-7-gb1e524c8c)
file : heat3d_small.edp
Load: lg_fem lg_mesh lg_mesh3 eigenvalue
(already loaded: msh3)
*** Warning The identifier y0 hide a Global identifier
sizestack + 1024 =4280 ( 3256 )

Time: 0.01 | min: -0.0282293 | max: 0.995666 | DMP violations: 1
Time: 0.02 | min: -0.0554627 | max: 0.986755 | DMP violations: 2
Time: 0.03 | min: -0.0814489 | max: 0.973381 | DMP violations: 3
Time: 0.04 | min: -0.106035 | max: 0.955699 | DMP violations: 4
Time: 0.05 | min: -0.129081 | max: 0.933897 | DMP violations: 5
Time: 0.06 | min: -0.150466 | max: 0.908196 | DMP violations: 6
Time: 0.07 | min: -0.170082 | max: 0.878846 | DMP violations: 7
Time: 0.08 | min: -0.18784 | max: 0.846127 | DMP violations: 8
Time: 0.09 | min: -0.20367 | max: 0.81034 | DMP violations: 9
Time: 0.1 | min: -0.217519 | max: 0.77181 | DMP violations: 10
Time: 0.11 | min: -0.229352 | max: 0.730877 | DMP violations: 11
Time: 0.12 | min: -0.239153 | max: 0.687897 | DMP violations: 12
Time: 0.13 | min: -0.246923 | max: 0.643235 | DMP violations: 13
Time: 0.14 | min: -0.252683 | max: 0.597262 | DMP violations: 14
Time: 0.15 | min: -0.256468 | max: 0.550353 | DMP violations: 15
Time: 0.16 | min: -0.258332 | max: 0.50288 | DMP violations: 16
Time: 0.17 | min: -0.258342 | max: 0.455211 | DMP violations: 17
Time: 0.18 | min: -0.256579 | max: 0.416945 | DMP violations: 18
Time: 0.19 | min: -0.253284 | max: 0.414025 | DMP violations: 19
Time: 0.2 | min: -0.251049 | max: 0.408155 | DMP violations: 20
Time: 0.21 | min: -0.24733 | max: 0.399426 | DMP violations: 21
Time: 0.22 | min: -0.242244 | max: 0.387955 | DMP violations: 22
Time: 0.23 | min: -0.235911 | max: 0.37389 | DMP violations: 23
Time: 0.24 | min: -0.22846 | max: 0.376792 | DMP violations: 24
Time: 0.25 | min: -0.22183 | max: 0.397626 | DMP violations: 25
Time: 0.26 | min: -0.226586 | max: 0.416604 | DMP violations: 26
Time: 0.27 | min: -0.230235 | max: 0.433504 | DMP violations: 27
Time: 0.28 | min: -0.232755 | max: 0.448126 | DMP violations: 28
Time: 0.29 | min: -0.234136 | max: 0.460299 | DMP violations: 29
Time: 0.3 | min: -0.234381 | max: 0.469879 | DMP violations: 30
Time: 0.31 | min: -0.233503 | max: 0.476752 | DMP violations: 31
Time: 0.32 | min: -0.231529 | max: 0.480835 | DMP violations: 32
Time: 0.33 | min: -0.228493 | max: 0.482077 | DMP violations: 33
Time: 0.34 | min: -0.224443 | max: 0.480462 | DMP violations: 34
Time: 0.35 | min: -0.219434 | max: 0.476002 | DMP violations: 35
Time: 0.36 | min: -0.213529 | max: 0.468744 | DMP violations: 36
Time: 0.37 | min: -0.206799 | max: 0.458764 | DMP violations: 37
Time: 0.38 | min: -0.199324 | max: 0.446169 | DMP violations: 38
Time: 0.39 | min: -0.191185 | max: 0.431095 | DMP violations: 39
Time: 0.4 | min: -0.18247 | max: 0.4137 | DMP violations: 40
Time: 0.41 | min: -0.17327 | max: 0.394171 | DMP violations: 41
Time: 0.42 | min: -0.163677 | max: 0.372711 | DMP violations: 42
Time: 0.43 | min: -0.153784 | max: 0.349544 | DMP violations: 43
Time: 0.44 | min: -0.143684 | max: 0.32491 | DMP violations: 44
Time: 0.45 | min: -0.13347 | max: 0.320886 | DMP violations: 45
Time: 0.46 | min: -0.142029 | max: 0.319538 | DMP violations: 46
Time: 0.47 | min: -0.151563 | max: 0.316236 | DMP violations: 47
Time: 0.48 | min: -0.159045 | max: 0.310992 | DMP violations: 48
Time: 0.49 | min: -0.164569 | max: 0.303836 | DMP violations: 49
Time: 0.5 | min: -0.168892 | max: 0.301431 | DMP violations: 50
Time: 0.51 | min: -0.171347 | max: 0.315554 | DMP violations: 51
Time: 0.52 | min: -0.17201 | max: 0.328413 | DMP violations: 52
Time: 0.53 | min: -0.170974 | max: 0.339825 | DMP violations: 53
Time: 0.54 | min: -0.168346 | max: 0.349616 | DMP violations: 54
Time: 0.55 | min: -0.168384 | max: 0.357622 | DMP violations: 55
Time: 0.56 | min: -0.176067 | max: 0.363692 | DMP violations: 56
Time: 0.57 | min: -0.182625 | max: 0.367691 | DMP violations: 57
Time: 0.58 | min: -0.188029 | max: 0.3695 | DMP violations: 58
Time: 0.59 | min: -0.192267 | max: 0.36902 | DMP violations: 59
Time: 0.6 | min: -0.195339 | max: 0.366174 | DMP violations: 60
Time: 0.61 | min: -0.197259 | max: 0.360904 | DMP violations: 61
Time: 0.62 | min: -0.198055 | max: 0.353177 | DMP violations: 62
Time: 0.63 | min: -0.197769 | max: 0.342985 | DMP violations: 63
Time: 0.64 | min: -0.196453 | max: 0.330344 | DMP violations: 64
Time: 0.65 | min: -0.194172 | max: 0.315292 | DMP violations: 65
Time: 0.66 | min: -0.190998 | max: 0.324328 | DMP violations: 66
Time: 0.67 | min: -0.187016 | max: 0.338559 | DMP violations: 67
Time: 0.68 | min: -0.182317 | max: 0.351959 | DMP violations: 68
Time: 0.69 | min: -0.176996 | max: 0.364402 | DMP violations: 69
Time: 0.7 | min: -0.171157 | max: 0.375771 | DMP violations: 70
Time: 0.71 | min: -0.164905 | max: 0.385953 | DMP violations: 71
Time: 0.72 | min: -0.158345 | max: 0.394843 | DMP violations: 72
Time: 0.73 | min: -0.151587 | max: 0.402342 | DMP violations: 73
Time: 0.74 | min: -0.144734 | max: 0.408365 | DMP violations: 74
Time: 0.75 | min: -0.137891 | max: 0.412835 | DMP violations: 75
Time: 0.76 | min: -0.14505 | max: 0.415687 | DMP violations: 76
Time: 0.77 | min: -0.162069 | max: 0.416869 | DMP violations: 77
Time: 0.78 | min: -0.178457 | max: 0.416343 | DMP violations: 78
Time: 0.79 | min: -0.194065 | max: 0.414083 | DMP violations: 79
Time: 0.8 | min: -0.208752 | max: 0.410079 | DMP violations: 80
Time: 0.81 | min: -0.222384 | max: 0.404336 | DMP violations: 81
Time: 0.82 | min: -0.234833 | max: 0.396872 | DMP violations: 82
Time: 0.83 | min: -0.245983 | max: 0.387721 | DMP violations: 83
Time: 0.84 | min: -0.255727 | max: 0.376931 | DMP violations: 84
Time: 0.85 | min: -0.263971 | max: 0.364565 | DMP violations: 85
Time: 0.86 | min: -0.270635 | max: 0.350866 | DMP violations: 86
Time: 0.87 | min: -0.292103 | max: 0.337605 | DMP violations: 87
Time: 0.88 | min: -0.311469 | max: 0.323062 | DMP violations: 88
Time: 0.89 | min: -0.328452 | max: 0.307337 | DMP violations: 89
Time: 0.9 | min: -0.342941 | max: 0.290538 | DMP violations: 90
Time: 0.91 | min: -0.354848 | max: 0.272785 | DMP violations: 91
Time: 0.92 | min: -0.364112 | max: 0.254201 | DMP violations: 92
Time: 0.93 | min: -0.370698 | max: 0.234919 | DMP violations: 93
Time: 0.94 | min: -0.374595 | max: 0.215075 | DMP violations: 94
Time: 0.95 | min: -0.37582 | max: 0.19481 | DMP violations: 95
Time: 0.96 | min: -0.374415 | max: 0.184427 | DMP violations: 96
Time: 0.97 | min: -0.370449 | max: 0.195187 | DMP violations: 97
Time: 0.98 | min: -0.364011 | max: 0.204447 | DMP violations: 98
Time: 0.99 | min: -0.355216 | max: 0.212092 | DMP violations: 99
Time: 1 | min: -0.344199 | max: 0.218027 | DMP violations: 100
Initial bounds: [2.66448e-147, 1]
Boundary value: 0
DMP bounds: [0, 1]
Total DMP violations: 100
FreeFem++.exe heat3d_small.edp -cd -wg -ns
– FreeFem++ v4.15 (Wed Dec 11 16:48:56 CET 2024 - git v4.15-7-gb1e524c8c)
file : heat3d_small.edp
Load: lg_fem lg_mesh lg_mesh3 eigenvalue
(already loaded: msh3)
*** Warning The identifier y0 hide a Global identifier
sizestack + 1024 =4280 ( 3256 )

Time: 0.01 | min: -0.0138431 | max: 0.999109 | DMP violations: 1
Time: 0.02 | min: -0.0277248 | max: 0.997268 | DMP violations: 2
Time: 0.03 | min: -0.0415371 | max: 0.994483 | DMP violations: 3
Time: 0.04 | min: -0.0552572 | max: 0.990759 | DMP violations: 4
Time: 0.05 | min: -0.0688626 | max: 0.986106 | DMP violations: 5
Time: 0.06 | min: -0.082331 | max: 0.980535 | DMP violations: 6
Time: 0.07 | min: -0.0956404 | max: 0.974059 | DMP violations: 7
Time: 0.08 | min: -0.108769 | max: 0.966694 | DMP violations: 8
Time: 0.09 | min: -0.121697 | max: 0.958457 | DMP violations: 9
Time: 0.1 | min: -0.134402 | max: 0.949366 | DMP violations: 10
Time: 0.11 | min: -0.146866 | max: 0.939443 | DMP violations: 11
Time: 0.12 | min: -0.159068 | max: 0.92871 | DMP violations: 12
Time: 0.13 | min: -0.170989 | max: 0.917192 | DMP violations: 13
Time: 0.14 | min: -0.182613 | max: 0.904914 | DMP violations: 14
Time: 0.15 | min: -0.193921 | max: 0.891905 | DMP violations: 15
Time: 0.16 | min: -0.204897 | max: 0.878193 | DMP violations: 16
Time: 0.17 | min: -0.215526 | max: 0.863809 | DMP violations: 17
Time: 0.18 | min: -0.225792 | max: 0.848785 | DMP violations: 18
Time: 0.19 | min: -0.235681 | max: 0.833153 | DMP violations: 19
Time: 0.2 | min: -0.245181 | max: 0.816947 | DMP violations: 20
Time: 0.21 | min: -0.25428 | max: 0.800204 | DMP violations: 21
Time: 0.22 | min: -0.262966 | max: 0.782957 | DMP violations: 22
Time: 0.23 | min: -0.271229 | max: 0.765246 | DMP violations: 23
Time: 0.24 | min: -0.279061 | max: 0.747106 | DMP violations: 24
Time: 0.25 | min: -0.286453 | max: 0.728575 | DMP violations: 25
Time: 0.26 | min: -0.293397 | max: 0.709693 | DMP violations: 26
Time: 0.27 | min: -0.299889 | max: 0.690498 | DMP violations: 27
Time: 0.28 | min: -0.305923 | max: 0.67103 | DMP violations: 28
Time: 0.29 | min: -0.311495 | max: 0.651327 | DMP violations: 29
Time: 0.3 | min: -0.316602 | max: 0.631428 | DMP violations: 30
Time: 0.31 | min: -0.321242 | max: 0.611374 | DMP violations: 31
Time: 0.32 | min: -0.325416 | max: 0.591202 | DMP violations: 32
Time: 0.33 | min: -0.329122 | max: 0.570952 | DMP violations: 33
Time: 0.34 | min: -0.332362 | max: 0.550662 | DMP violations: 34
Time: 0.35 | min: -0.335139 | max: 0.53037 | DMP violations: 35
Time: 0.36 | min: -0.337456 | max: 0.510113 | DMP violations: 36
Time: 0.37 | min: -0.339316 | max: 0.489928 | DMP violations: 37
Time: 0.38 | min: -0.340726 | max: 0.46985 | DMP violations: 38
Time: 0.39 | min: -0.341691 | max: 0.449914 | DMP violations: 39
Time: 0.4 | min: -0.342218 | max: 0.430155 | DMP violations: 40
Time: 0.41 | min: -0.342315 | max: 0.410604 | DMP violations: 41
Time: 0.42 | min: -0.341991 | max: 0.391294 | DMP violations: 42
Time: 0.43 | min: -0.341255 | max: 0.372256 | DMP violations: 43
Time: 0.44 | min: -0.340117 | max: 0.353518 | DMP violations: 44
Time: 0.45 | min: -0.338588 | max: 0.335108 | DMP violations: 45
Time: 0.46 | min: -0.336681 | max: 0.317054 | DMP violations: 46
Time: 0.47 | min: -0.334406 | max: 0.29938 | DMP violations: 47
Time: 0.48 | min: -0.331777 | max: 0.28211 | DMP violations: 48
Time: 0.49 | min: -0.328807 | max: 0.276079 | DMP violations: 49
Time: 0.5 | min: -0.32551 | max: 0.277552 | DMP violations: 50
Time: 0.51 | min: -0.3219 | max: 0.278577 | DMP violations: 51
Time: 0.52 | min: -0.317993 | max: 0.279153 | DMP violations: 52
Time: 0.53 | min: -0.313802 | max: 0.279281 | DMP violations: 53
Time: 0.54 | min: -0.309345 | max: 0.278962 | DMP violations: 54
Time: 0.55 | min: -0.304635 | max: 0.2782 | DMP violations: 55
Time: 0.56 | min: -0.299689 | max: 0.277 | DMP violations: 56
Time: 0.57 | min: -0.295305 | max: 0.275368 | DMP violations: 57
Time: 0.58 | min: -0.297038 | max: 0.273313 | DMP violations: 58
Time: 0.59 | min: -0.298458 | max: 0.270843 | DMP violations: 59
Time: 0.6 | min: -0.299566 | max: 0.267968 | DMP violations: 60
Time: 0.61 | min: -0.300364 | max: 0.264702 | DMP violations: 61
Time: 0.62 | min: -0.300854 | max: 0.261056 | DMP violations: 62
Time: 0.63 | min: -0.301039 | max: 0.257045 | DMP violations: 63
Time: 0.64 | min: -0.300923 | max: 0.252684 | DMP violations: 64
Time: 0.65 | min: -0.300512 | max: 0.24799 | DMP violations: 65
Time: 0.66 | min: -0.29981 | max: 0.242979 | DMP violations: 66
Time: 0.67 | min: -0.298823 | max: 0.237671 | DMP violations: 67
Time: 0.68 | min: -0.29756 | max: 0.232083 | DMP violations: 68
Time: 0.69 | min: -0.296026 | max: 0.226235 | DMP violations: 69
Time: 0.7 | min: -0.294231 | max: 0.220148 | DMP violations: 70
Time: 0.71 | min: -0.292182 | max: 0.213842 | DMP violations: 71
Time: 0.72 | min: -0.289889 | max: 0.207339 | DMP violations: 72
Time: 0.73 | min: -0.287362 | max: 0.20066 | DMP violations: 73
Time: 0.74 | min: -0.28461 | max: 0.193828 | DMP violations: 74
Time: 0.75 | min: -0.281644 | max: 0.186864 | DMP violations: 75
Time: 0.76 | min: -0.278475 | max: 0.17979 | DMP violations: 76
Time: 0.77 | min: -0.275113 | max: 0.172629 | DMP violations: 77
Time: 0.78 | min: -0.271572 | max: 0.165403 | DMP violations: 78
Time: 0.79 | min: -0.267861 | max: 0.158133 | DMP violations: 79
Time: 0.8 | min: -0.263994 | max: 0.150843 | DMP violations: 80
Time: 0.81 | min: -0.259982 | max: 0.147974 | DMP violations: 81
Time: 0.82 | min: -0.255837 | max: 0.147477 | DMP violations: 82
Time: 0.83 | min: -0.251571 | max: 0.146701 | DMP violations: 83
Time: 0.84 | min: -0.247198 | max: 0.145643 | DMP violations: 84
Time: 0.85 | min: -0.242729 | max: 0.147074 | DMP violations: 85
Time: 0.86 | min: -0.238176 | max: 0.148963 | DMP violations: 86
Time: 0.87 | min: -0.233552 | max: 0.150707 | DMP violations: 87
Time: 0.88 | min: -0.228869 | max: 0.152301 | DMP violations: 88
Time: 0.89 | min: -0.224139 | max: 0.153743 | DMP violations: 89
Time: 0.9 | min: -0.219374 | max: 0.155031 | DMP violations: 90
Time: 0.91 | min: -0.214585 | max: 0.156162 | DMP violations: 91
Time: 0.92 | min: -0.209783 | max: 0.157134 | DMP violations: 92
Time: 0.93 | min: -0.20498 | max: 0.157947 | DMP violations: 93
Time: 0.94 | min: -0.201545 | max: 0.1586 | DMP violations: 94
Time: 0.95 | min: -0.198245 | max: 0.159093 | DMP violations: 95
Time: 0.96 | min: -0.194783 | max: 0.159426 | DMP violations: 96
Time: 0.97 | min: -0.191167 | max: 0.159599 | DMP violations: 97
Time: 0.98 | min: -0.187763 | max: 0.159615 | DMP violations: 98
Time: 0.99 | min: -0.189905 | max: 0.159474 | DMP violations: 99
Time: 1 | min: -0.19176 | max: 0.159179 | DMP violations: 100
Initial bounds: [2.66448e-147, 1]
Boundary value: 0
DMP bounds: [0, 1]
Total DMP violations: 100

2

Your code gives the discrete maximum principle satisfied with mass lumping if you don’t put advection (bx=by=bz=0). It you want it to work also with advection, I think you need some kind of upwinding.

1 Like
3

Even when I use convect, it doesn’t work. I tried in 2D; convection gives a positive solution. But in 3D, it still violates the DMP. I suspect it’s the mesh. Do you think this is true?

4

I think there is a problem with convect in 3d.
With the following test extracted from your code

verbosity = 0;

// Mesh and time parameters
int nx = 8;
int ny = 8;
int nz = 8;
real dt = 0.001;
real bx = 1.0;
real by = 0.5;
real bz = 0.5;


mesh3 Th = cube(nx, ny, nz);
fespace Vh(Th, P1);
Vh uold, btrx,btry,btrz;
btrx=bx;
btry=by;
btrz=bz;
real xx0 = 0.25;
real yy0 = 0.25;
real zz0 = 0.25;
real sigma = 0.05;
uold = exp(-((x - xx0)^2 + (y - yy0)^2 + (z - zz0)^2) / (2.0 * sigma^2));

Vh ueval;
ueval=convect([btrx,btry,btrz],-dt,uold);
cout << "uoldmin " << uold[].min << " uevalmin " << ueval[].min << endl;
for (int ii=0;ii<Vh.ndof;ii++){
 if (ueval[][ii]<0.) {
  cout << "position [" << Th(ii).x << ", " << Th(ii).y << ", " << Th(ii).z << "]" ;
  cout << " ueval=" << ueval[][ii] << endl;
 }
}

the convection of a positive function gives negative values.

1 Like
5

Yes. Do you have any suggestions on how to solve this issue?

Thank you.