verbosity = 0; // to limit message printing load "Element_P3"; // load 3rd-order elements, used for velocity load "Element_P4"; // load 4th-order elements, used for stream function //mesh resolution - mesh points per unit of boundary - no need to change int n=6; //parameters - NEED TO ADD GRAD-DIV PARAMETER real nu = 1.0e-2; // kinematic viscosity // penalty method for uniqueness of pressure; needed for the linear solve real penalty = 1.0e-11; // I do not recommend changing this. //set parameters for nonlinear iteration - tread carefully if you change these int nliter, maxiter=300; real nlerr, nltol = 1.0e-9; // manufactured solutions - DO NOT CHANGE! real amp=1.0/8000; // scales solution so terms of NSE are not too large func u1True=amp*(2*y-1+(20-x)^3); func u2True=amp*(3*y*(20-x)^2 -2*x); func pTrue=amp*(20-x)*(2*y-1); //RHS functions that force the manufactured solutions - DO NOT CHANGE! func f1=amp*(2*amp*(3*y*(20-x)^2-2*x)+1-2*y-6*nu*(20-x)-3*amp*((20-x)^2)*(2*y-1+(20-x)^3)); func f2=amp*(40-2*x-2*amp*((20-x)^3+2*y-1)*(1+3*y*(20-x))-6*y*nu+3*amp*((20-x)^2)*(3*y*(20-x)^2-2*x)); //domain specifications - DO NOT CHANGE! mesh Omega; // **************************** C3 ********************************* // * * // C4 * // * * // *******C5******* Omega C2 // * * // C6 * // * * // ************* C1 ********************************* border C1(t=5,20){ x=t; y=0; label=1; } border C2(t=0,1){ x=20; y=t; label=2; } border C3(t=0,20) { x=20-t; y=1; label=1; } border C4(t=0,1) { x=0; y=1-0.5*t; label=3; } border C5(t=0,5) { x=t; y=0.5; label=1; } border C6(t=0,1) { x=5; y=0.5-0.5*t; label=1; } // label = 1 means the top and bottom walls // label = 2 means the outflow // label = 3 means the inflow Omega=buildmesh(C1(15*n) + C2(n) + C3(20*n) + C4(n/2)+C5(5*n)+C6(n/2)); //savemesh(Omega,"Omega_"+n+".msh"); // Uncomment the following to see the mesh. // It will wait for you to hit enter/return before solving the problem. // plot(Omega,wait=1); //finite element spaces - DO NOT CHANGE! fespace X(Omega,P3); // cubics for velocity components fespace Q(Omega,P2); // quads for pressure fespace S(Omega,P4); // for stream function int nn = X.ndof, nnn = Q.ndof; // degrees of freedom for each space cout << "Degrees of freedom over all subdomains" << endl; cout << "Velocity: " << 2*nn << " Pressure: " << nnn << endl; // Define variables and test functions - NO NEED TO CHANGE OR ADD HERE X u1n,u2n,u1nplus1,u2nplus1,v1,v2; Q pn,q1; S psi,phi; // diffusion terms - GOOD PLACE TO ADD GRAD-DIV TERMS varf a11a(u1nplus1,v1) = int2d(Omega)(nu*(2*dx(u1nplus1)*dx(v1)+dy(u1nplus1)*dy(v1))) +on(1,3,u1nplus1=u1True); varf a12a(u1nplus1,v2) = int2d(Omega)(nu*dy(u1nplus1)*dx(v2)); varf a21a(u2nplus1,v1) = int2d(Omega)(nu*dx(u2nplus1)*dy(v1)); varf a22a(u2nplus1,v2) = int2d(Omega)(nu*(dx(u2nplus1)*dx(v2)+2*dy(u2nplus1)*dy(v2))) +on(1,3,u2nplus1=u2True); // divergence constraint - DO NOT CHANGE! varf div1(u1nplus1,q1) = int2d(Omega)(1.0*dx(u1nplus1)*q1); varf div2(u2nplus1,q1) = int2d(Omega)(1.0*dy(u2nplus1)*q1); // pressure terms - DO NOT CHANGE! varf pres1(pn,v1) = int2d(Omega)(-1.0*pn*dx(v1)); varf pres2(pn,v2) = int2d(Omega)(-1.0*pn*dy(v2)); varf pq1(pn,q1) = int2d(Omega)(-1.0*penalty*pn*q1); //right-hand side forms varf ff1(u1nplus1,v1) = int2d(Omega,qft=qf9pT)(f1*v1) +on(1,3,u1nplus1=u1True); varf ff2(u2nplus1,v2) = int2d(Omega,qft=qf9pT)(f2*v2) +on(1,3,u2nplus1=u2True); varf gg1(pn,q1) = int2d(Omega)(0.0*q1); // divergence = 0 //Variational forms that are dependent on nonlinear iteration /* ADD NONLINEARITY TERMS HERE */ //allocate the solution and right-hand side vectors - DO NOT CHANGE! real[int] F1(nn), F2(nn), F(2*nn+nnn), G1(nnn), Usol(2*nn+nnn), Uoldsol(2*nn+nnn), Usoldiff(2*nn+nnn); /* initial values - you should get the right answer with 1 iteration if you use the true solution here */ u1n = u1True; u2n = u2True; pn = pTrue; /* u1n = 0; u2n = 0; pn = 0; */ //assemble matrix blocks that are indepenent of the nonlinear iteration /* No changes needed here unless you defined new "varf" structures for the grad-div terms */ matrix A11a = a11a(X,X); matrix A12a = a12a(X,X); matrix A21a = a21a(X,X); matrix A22a = a22a(X,X); matrix Div1 = div1(X,Q); matrix Div2 = div2(X,Q); matrix Pres1 = pres1(Q,X); matrix Pres2 = pres2(Q,X); matrix PQ1 = pq1(Q,Q); G1 = gg1(0,Q); F1 = ff1(0,X); F2 = ff2(0,X); //compute the rhs vectors - should not need to do anything here F = [F1,F2,G1]; //begin nonlinear iteration for(nliter=1;nliter<=maxiter;nliter++) { Uoldsol = [u1n[], u2n[], pn[]]; // ADD IN NONLINEAR MATRIX CONTRIBUTIONS, FORM A11b, A22b HERE! //assemble the coefficient matrix - DO NOT CHANGE! matrix A = [[A11b, A21a, Pres1], [A12a, A22b, Pres2], [Div1, Div2, PQ1]]; //solve the problem - DO NOT CHANGE! set(A, solver=UMFPACK); Usol = A^-1*F; //get the solution components - DO NOT CHANGE! [u1nplus1[],u2nplus1[],pn[]] = Usol; //check for convergence - DO NOT CHANGE! Usoldiff = Usol-Uoldsol; nlerr = Usoldiff.l2/Usol.l2; cout << "Nonlinear iteration error: " << nlerr << endl; if (nlerr < nltol || nlerr > 1.0e6) break; - DO NOT CHANGE! //update temp solution - DO NOT CHANGE! u1n = u1nplus1; u2n = u2nplus1; Uoldsol = Usol; // plot solution each iteration if you want, otherwise comment it out plot(pn,[u1n,u2n],wait=0,fill=0,value=1,cmm="Computed Approximation"); }//end nonlinear iteration cout << "Nonlinear iterations required: " << nliter << endl; // capture error and output to screen - DO NOT CHANGE! real u1L2 = sqrt(int2d(Omega,qft=qf7pT)((u1n-u1True)^2)); cout<<"Horizontal velocity error = "<