diff --git a/Solver/CMakeLists.txt b/Solver/CMakeLists.txt
index a7701284f5148266b0e3e1c1a97360c381d8bda7..4bff04102c84113b4469be7f5960d5bc488451dd 100644
--- a/Solver/CMakeLists.txt
+++ b/Solver/CMakeLists.txt
@@ -19,7 +19,7 @@ set(SRC
dgConservationLawShallowWater2d.cpp
dgConservationLawWaveEquation.cpp
dgConservationLawPerfectGas.cpp
- SlopeLimiter.cpp
+ dgSlopeLimiter.cpp
function.cpp
functionDerivator.cpp
luaFunction.cpp
diff --git a/Solver/SlopeLimiter.cpp b/Solver/SlopeLimiter.cpp
deleted file mode 100644
index f8f3f56cb7be073b06fd8b2db4a6c6a33b63607c..0000000000000000000000000000000000000000
--- a/Solver/SlopeLimiter.cpp
+++ /dev/null
@@ -1,139 +0,0 @@
-#include "dgLimiter.h"
-#include "dgGroupOfElements.h"
-#include "dgSystemOfEquations.h"
-
-
-// //----------------------------------------------------------------------------------
-// bool SlopeLimiter::apply ( dgDofContainer &solution,
-// std::vector<dgGroupOfElements*> &eGroups,
-// std::vector<dgGroupOfFaces*> &fGroups)
-// {
-
-// //WARNING FOR ONLY 1 GROUP OF FACES
-// dgGroupOfFaces* group = fGroups[0];
-// fullMatrix<double> SolLeft = solution._dataProxys[0];
-// fullMatrix<double> SolRight = solution._dataProxys[0];
-// int nbFields = solution->getNbFields();
-// int totNbElems = solution->getNbElements();
-
-// // first compute max and min of all fields for all stencils
-// //----------------------------------------------------------
-// fullMatrix<double> MIN (totNbElems,nbFields );
-// fullMatrix<double> MAX (totNbElems, nbFields);
-// MIN.set_all ( 1.e22);
-// MAX.set_all (-1.e22);
-
-// fullMatrix<double> MEANL (totNbElems, nbFields);
-// fullMatrix<double> MEANR (totNbElems, nbFields);
-// MEANL.set_all ( 0.01);
-// MEANR.set_all ( 0.01);
-
-// for(int iFace=0 ; iFace<group->getNbElements();iFace++)
-// {
-// int iElementL = group->getElementLeftId(iFace);
-// int iElementR = group->getElementRightId(iFace);
-// if (iElementR >= 0)
-// {
-// fullMatrix<double> TempL, TempR;
-// TempL.setAsProxy(SolLeft, nbFields*iELementL, nbFields );
-// TempR.setAsProxy(SolRight, nbFields*iELementR, nbFields );
-// printf("TempL size 1 = %d %d ", TempL.size1(), TempL.size2());
-// exit(1);
-
-// int fSize = TempL.size2();
-// for (int k=0; k< fSize; ++k)
-// {
-// double AVGL = 0;
-// double AVGR = 0;
-// for (int i=0; i<fSize; ++i)
-// {
-// AVGL += TempL(i,k);
-// AVGR += TempR(i,k);
-// }
-// AVGL /= (double) fSize;
-// AVGR /= (double) fSize;
-// MIN ( iElementL , k ) = std::min ( AVGR , MIN ( iElementL , k ) );
-// MAX ( iElementL , k ) = std::max ( AVGR , MAX ( iElementL , k ) );
-// MIN ( iElementR , k ) = std::min ( AVGL , MIN ( iElementR , k ) );
-// MAX ( iElementR , k ) = std::max ( AVGL , MAX ( iElementR , k ) );
-
-// MEANR( iElementL,k ) = AVGR;
-// MEANL( iElementR,k ) = AVGL;
-
-// }
-// }
-// else{
-// fullMatrix<double> TempL, ;
-// TempL.setAsProxy(SolLeft, nbFields*iELementL, nbFields );
-// for (int k=0; k<nbFields; ++k){
-// for (int i=0; i<fSize; ++i){
-// MIN ( iElementL , k ) = std::min ( TempL(i,k) , MIN ( iElementL , k ) );
-// MAX ( iElementL , k ) = std::max ( TempL(i,k) , MAX ( iElementL , k ) );
-
-
-// double AVG = 0;
-// for (int i=0; i<fSize; ++i) AVG += TempL(i,k);
-
-// }
-// }
-// }
-// }
-// // some parallel should be done here in order to compute averages on other processors
-// //----------------------------------------------------------
-
-// //...
-
-// // then limit the solution
-// //----------------------------------------------------------
-
-// // for (int iElement=0 ; iElement<totNbElems; ++iElement)
-// // {
-// // fullMatrix Temp = fullMatrix(solution->touchElementCoefs (iElement));
-// // for (int k=0; k<nbFields; ++k)
-// // {
-// // double AVG = 0.;
-// // double locMax = -1.e22;
-// // double locMin = 1.e22;
-// // double neighMax = MAX (iElement,k);
-// // double neighMin = MIN (iElement,k);
-// // for (int i=0; i<fSize; ++i)
-// // {
-// // AVG += Temp(i,k);
-// // locMax = std::max (locMax, Temp (i,k));
-// // locMin = std::min (locMin, Temp (i,k));
-// // }
-// // AVG /= (double) fSize;
-
-// // //SLOPE LIMITING DG
-// // //-------------------
-// // for (int i=0; i<fSize; ++i)Temp(i,k) -= AVG;
-
-// // double slopeLimiterValue = 1.0;
-// // if (locMax != AVG && locMax > neighMax) slopeLimiterValue = (neighMax-AVG) / (locMax-AVG);
-// // if (locMin != AVG && locMin < neighMin) slopeLimiterValue = std::min ( slopeLimiterValue , (AVG-neighMin) / (AVG-locMin) );
-// // if (AVG < neighMin) slopeLimiterValue = 0;
-// // if (AVG > neighMax) slopeLimiterValue = 0;
-
-// // if (detectDISC(iElement) == 0.0) slopeLimiterValue = 1.0; // do not limit
-// // //if (slopeLimiterValue != 1.0) printf("limit (%g) elem =%d \n", slopeLimiterValue, iElement);
-
-// // for (int i=0; i<fSize; ++i) Temp(i,k) *= slopeLimiterValue;
-
-// // for (int i=0; i<fSize; ++i) Temp(i,k) += AVG;
-
-// // //MINMOD DG (IF CHANGE TO THIS ADD LIMITER LINE IN RUNGEKUTTANNOPERATOR.cc)
-// // //--------------------------------
-// // // if (detectDISC(iElement) != 0.0){
-// // // const size_t nn = Temp.size1();
-// // // double alpha = 0.5; //0.5 1.0 1.3 (alpha=0.5 => MUSC Schem VAN LEER)
-// // // Temp(0,k) = AVG-minmod(AVG-Temp(0,k), alpha*(AVG-MEANL(iElement,k)), alpha*(MEANR(iElement,k)-AVG));
-// // // Temp(nn-1,k) = AVG+minmod(Temp(nn-1,k)-AVG, alpha*(AVG-MEANL(iElement,k)), alpha*(MEANR(iElement,k)-AVG));
-// // // }
-
-
-// // }
-// // }
-
-// return true;
-
-// }
diff --git a/Solver/TESTCASES/Advection1D.lua b/Solver/TESTCASES/Advection1D.lua
index f0fe33c73a4457fc495a7c6fb5827ffa627f03aa..fb2abf5186f4023214c73a0e56b7e0d634e0b75e 100644
--- a/Solver/TESTCASES/Advection1D.lua
+++ b/Solver/TESTCASES/Advection1D.lua
@@ -1,7 +1,7 @@
model = GModel ()
model:load ('edge.geo')
---model:mesh (1)
---model:save ('edge_new.msh')
+-- model:mesh (1)
+-- model:save ('edge.msh')
model:load ('edge.msh')
dg = dgSystemOfEquations (model)
dg:setOrder(1)
@@ -37,17 +37,24 @@ function initial_condition( xyz , f )
x = xyz:get(i,0)
y = xyz:get(i,1)
z = xyz:get(i,2)
- f:set (i, 0, math.exp(-100*((x+0.8)^2)))
+ if (x>-0.3 and x<0.3) then
+ f:set (i, 0, 1)
+ else
+ f:set (i, 0, 0)
+ end
end
end
dg:L2Projection(FunctionLua(1,'initial_condition',{'XYZ'}):getName())
+dg:exportSolution('output/Adv1D_unlimited')
+dg:limitSolution()
dg:exportSolution('output/Adv1D_00000')
+
-- main loop
n = 5
for i=1,100*n do
- norm = dg:RK44(0.03)
+ norm = dg:RK44_limiter(0.03)
if (i % n == 0) then
print('iter',i,norm)
dg:exportSolution(string.format("output/Adv1D-%05d", i))
diff --git a/Solver/TESTCASES/edge.msh b/Solver/TESTCASES/edge.msh
new file mode 100644
index 0000000000000000000000000000000000000000..5ee49d3930ffa1973e9913f8c37f65abe2d39865
--- /dev/null
+++ b/Solver/TESTCASES/edge.msh
@@ -0,0 +1,416 @@
+$MeshFormat
+2.1 0 8
+$EndMeshFormat
+$PhysicalNames
+3
+0 1 "Left"
+0 2 "Right"
+1 3 "Line"
+$EndPhysicalNames
+$Nodes
+200
+1 -1 0 0
+2 1 0 0
+3 -0.9899497487437465 0 0
+4 -0.9798994974874929 0 0
+5 -0.9698492462312395 0 0
+6 -0.959798994974986 0 0
+7 -0.9497487437187324 0 0
+8 -0.9396984924624789 0 0
+9 -0.9296482412062255 0 0
+10 -0.9195979899499719 0 0
+11 -0.9095477386937184 0 0
+12 -0.8994974874374648 0 0
+13 -0.8894472361812114 0 0
+14 -0.8793969849249579 0 0
+15 -0.8693467336687044 0 0
+16 -0.8592964824124508 0 0
+17 -0.8492462311561974 0 0
+18 -0.8391959798999438 0 0
+19 -0.8291457286436903 0 0
+20 -0.8190954773874368 0 0
+21 -0.8090452261311832 0 0
+22 -0.7989949748749298 0 0
+23 -0.7889447236186762 0 0
+24 -0.7788944723624227 0 0
+25 -0.7688442211061692 0 0
+26 -0.7587939698499158 0 0
+27 -0.7487437185936622 0 0
+28 -0.7386934673374087 0 0
+29 -0.7286432160811551 0 0
+30 -0.7185929648249016 0 0
+31 -0.7085427135686482 0 0
+32 -0.6984924623123947 0 0
+33 -0.6884422110561411 0 0
+34 -0.6783919597998875 0 0
+35 -0.6683417085436341 0 0
+36 -0.6582914572873806 0 0
+37 -0.6482412060311271 0 0
+38 -0.6381909547748736 0 0
+39 -0.6281407035186201 0 0
+40 -0.6180904522623665 0 0
+41 -0.608040201006113 0 0
+42 -0.5979899497498595 0 0
+43 -0.5879396984936061 0 0
+44 -0.5778894472373525 0 0
+45 -0.567839195981099 0 0
+46 -0.5577889447248454 0 0
+47 -0.5477386934685919 0 0
+48 -0.5376884422123385 0 0
+49 -0.527638190956085 0 0
+50 -0.5175879396998314 0 0
+51 -0.5075376884435778 0 0
+52 -0.4974874371873209 0 0
+53 -0.4874371859310535 0 0
+54 -0.477386934674786 0 0
+55 -0.4673366834185185 0 0
+56 -0.4572864321622511 0 0
+57 -0.4472361809059836 0 0
+58 -0.4371859296497161 0 0
+59 -0.4271356783934487 0 0
+60 -0.4170854271371812 0 0
+61 -0.4070351758809138 0 0
+62 -0.3969849246246463 0 0
+63 -0.3869346733683788 0 0
+64 -0.3768844221121114 0 0
+65 -0.3668341708558439 0 0
+66 -0.3567839195995764 0 0
+67 -0.346733668343309 0 0
+68 -0.3366834170870415 0 0
+69 -0.3266331658307741 0 0
+70 -0.3165829145745066 0 0
+71 -0.3065326633182391 0 0
+72 -0.2964824120619717 0 0
+73 -0.2864321608057042 0 0
+74 -0.2763819095494368 0 0
+75 -0.2663316582931693 0 0
+76 -0.2562814070369018 0 0
+77 -0.2462311557806344 0 0
+78 -0.2361809045243669 0 0
+79 -0.2261306532680994 0 0
+80 -0.216080402011832 0 0
+81 -0.2060301507555645 0 0
+82 -0.1959798994992971 0 0
+83 -0.1859296482430296 0 0
+84 -0.1758793969867621 0 0
+85 -0.1658291457304947 0 0
+86 -0.1557788944742272 0 0
+87 -0.1457286432179598 0 0
+88 -0.1356783919616923 0 0
+89 -0.1256281407054248 0 0
+90 -0.1155778894491574 0 0
+91 -0.1055276381928899 0 0
+92 -0.09547738693662244 0 0
+93 -0.08542713568035498 0 0
+94 -0.07537688442408752 0 0
+95 -0.06532663316782006 0 0
+96 -0.0552763819115526 0 0
+97 -0.04522613065528525 0 0
+98 -0.03517587939901778 0 0
+99 -0.02512562814275032 0 0
+100 -0.01507537688648286 0 0
+101 -0.005025125630215399 0 0
+102 0.005025125626066052 0 0
+103 0.01507537688236149 0 0
+104 0.02512562813865671 0 0
+105 0.03517587939495215 0 0
+106 0.04522613065124759 0 0
+107 0.0552763819075428 0 0
+108 0.06532663316383824 0 0
+109 0.07537688442013346 0 0
+110 0.0854271356764289 0 0
+111 0.09547738693272434 0 0
+112 0.1055276381890196 0 0
+113 0.115577889445315 0 0
+114 0.1256281407016102 0 0
+115 0.1356783919579057 0 0
+116 0.1457286432142011 0 0
+117 0.1557788944704963 0 0
+118 0.1658291457267917 0 0
+119 0.1758793969830872 0 0
+120 0.1859296482393824 0 0
+121 0.1959798994956778 0 0
+122 0.2060301507519731 0 0
+123 0.2160804020082685 0 0
+124 0.2261306532645639 0 0
+125 0.2361809045208592 0 0
+126 0.2462311557771546 0 0
+127 0.25628140703345 0 0
+128 0.2663316582897453 0 0
+129 0.2763819095460407 0 0
+130 0.2864321608023359 0 0
+131 0.2964824120586314 0 0
+132 0.3065326633149268 0 0
+133 0.316582914571222 0 0
+134 0.3266331658275174 0 0
+135 0.3366834170838127 0 0
+136 0.3467336683401081 0 0
+137 0.3567839195964035 0 0
+138 0.3668341708526988 0 0
+139 0.3768844221089942 0 0
+140 0.3869346733652896 0 0
+141 0.3969849246215849 0 0
+142 0.4070351758778803 0 0
+143 0.4170854271341755 0 0
+144 0.427135678390471 0 0
+145 0.4371859296467664 0 0
+146 0.4472361809030616 0 0
+147 0.457286432159357 0 0
+148 0.4673366834156525 0 0
+149 0.4773869346719477 0 0
+150 0.4874371859282431 0 0
+151 0.4974874371845384 0 0
+152 0.5075376884408442 0 0
+153 0.5175879396971537 0 0
+154 0.5276381909534629 0 0
+155 0.5376884422097721 0 0
+156 0.5477386934660815 0 0
+157 0.5577889447223907 0 0
+158 0.5678391959787001 0 0
+159 0.5778894472350093 0 0
+160 0.5879396984913188 0 0
+161 0.597989949747628 0 0
+162 0.6080402010039372 0 0
+163 0.6180904522602466 0 0
+164 0.6281407035165558 0 0
+165 0.6381909547728652 0 0
+166 0.6482412060291745 0 0
+167 0.6582914572854837 0 0
+168 0.6683417085417931 0 0
+169 0.6783919597981023 0 0
+170 0.6884422110544117 0 0
+171 0.6984924623107209 0 0
+172 0.7085427135670304 0 0
+173 0.7185929648233396 0 0
+174 0.7286432160796488 0 0
+175 0.7386934673359582 0 0
+176 0.7487437185922674 0 0
+177 0.7587939698485768 0 0
+178 0.768844221104886 0 0
+179 0.7788944723611955 0 0
+180 0.7889447236175047 0 0
+181 0.7989949748738139 0 0
+182 0.8090452261301233 0 0
+183 0.8190954773864325 0 0
+184 0.8291457286427419 0 0
+185 0.8391959798990511 0 0
+186 0.8492462311553606 0 0
+187 0.8592964824116698 0 0
+188 0.869346733667979 0 0
+189 0.8793969849242884 0 0
+190 0.8894472361805976 0 0
+191 0.8994974874369071 0 0
+192 0.9095477386932163 0 0
+193 0.9195979899495255 0 0
+194 0.9296482412058349 0 0
+195 0.9396984924621441 0 0
+196 0.9497487437184535 0 0
+197 0.9597989949747627 0 0
+198 0.9698492462310722 0 0
+199 0.9798994974873814 0 0
+200 0.9899497487436906 0 0
+$EndNodes
+$Elements
+201
+1 15 3 1 1 0 1
+2 15 3 2 2 0 2
+3 1 3 3 1 0 1 3
+4 1 3 3 1 0 3 4
+5 1 3 3 1 0 4 5
+6 1 3 3 1 0 5 6
+7 1 3 3 1 0 6 7
+8 1 3 3 1 0 7 8
+9 1 3 3 1 0 8 9
+10 1 3 3 1 0 9 10
+11 1 3 3 1 0 10 11
+12 1 3 3 1 0 11 12
+13 1 3 3 1 0 12 13
+14 1 3 3 1 0 13 14
+15 1 3 3 1 0 14 15
+16 1 3 3 1 0 15 16
+17 1 3 3 1 0 16 17
+18 1 3 3 1 0 17 18
+19 1 3 3 1 0 18 19
+20 1 3 3 1 0 19 20
+21 1 3 3 1 0 20 21
+22 1 3 3 1 0 21 22
+23 1 3 3 1 0 22 23
+24 1 3 3 1 0 23 24
+25 1 3 3 1 0 24 25
+26 1 3 3 1 0 25 26
+27 1 3 3 1 0 26 27
+28 1 3 3 1 0 27 28
+29 1 3 3 1 0 28 29
+30 1 3 3 1 0 29 30
+31 1 3 3 1 0 30 31
+32 1 3 3 1 0 31 32
+33 1 3 3 1 0 32 33
+34 1 3 3 1 0 33 34
+35 1 3 3 1 0 34 35
+36 1 3 3 1 0 35 36
+37 1 3 3 1 0 36 37
+38 1 3 3 1 0 37 38
+39 1 3 3 1 0 38 39
+40 1 3 3 1 0 39 40
+41 1 3 3 1 0 40 41
+42 1 3 3 1 0 41 42
+43 1 3 3 1 0 42 43
+44 1 3 3 1 0 43 44
+45 1 3 3 1 0 44 45
+46 1 3 3 1 0 45 46
+47 1 3 3 1 0 46 47
+48 1 3 3 1 0 47 48
+49 1 3 3 1 0 48 49
+50 1 3 3 1 0 49 50
+51 1 3 3 1 0 50 51
+52 1 3 3 1 0 51 52
+53 1 3 3 1 0 52 53
+54 1 3 3 1 0 53 54
+55 1 3 3 1 0 54 55
+56 1 3 3 1 0 55 56
+57 1 3 3 1 0 56 57
+58 1 3 3 1 0 57 58
+59 1 3 3 1 0 58 59
+60 1 3 3 1 0 59 60
+61 1 3 3 1 0 60 61
+62 1 3 3 1 0 61 62
+63 1 3 3 1 0 62 63
+64 1 3 3 1 0 63 64
+65 1 3 3 1 0 64 65
+66 1 3 3 1 0 65 66
+67 1 3 3 1 0 66 67
+68 1 3 3 1 0 67 68
+69 1 3 3 1 0 68 69
+70 1 3 3 1 0 69 70
+71 1 3 3 1 0 70 71
+72 1 3 3 1 0 71 72
+73 1 3 3 1 0 72 73
+74 1 3 3 1 0 73 74
+75 1 3 3 1 0 74 75
+76 1 3 3 1 0 75 76
+77 1 3 3 1 0 76 77
+78 1 3 3 1 0 77 78
+79 1 3 3 1 0 78 79
+80 1 3 3 1 0 79 80
+81 1 3 3 1 0 80 81
+82 1 3 3 1 0 81 82
+83 1 3 3 1 0 82 83
+84 1 3 3 1 0 83 84
+85 1 3 3 1 0 84 85
+86 1 3 3 1 0 85 86
+87 1 3 3 1 0 86 87
+88 1 3 3 1 0 87 88
+89 1 3 3 1 0 88 89
+90 1 3 3 1 0 89 90
+91 1 3 3 1 0 90 91
+92 1 3 3 1 0 91 92
+93 1 3 3 1 0 92 93
+94 1 3 3 1 0 93 94
+95 1 3 3 1 0 94 95
+96 1 3 3 1 0 95 96
+97 1 3 3 1 0 96 97
+98 1 3 3 1 0 97 98
+99 1 3 3 1 0 98 99
+100 1 3 3 1 0 99 100
+101 1 3 3 1 0 100 101
+102 1 3 3 1 0 101 102
+103 1 3 3 1 0 102 103
+104 1 3 3 1 0 103 104
+105 1 3 3 1 0 104 105
+106 1 3 3 1 0 105 106
+107 1 3 3 1 0 106 107
+108 1 3 3 1 0 107 108
+109 1 3 3 1 0 108 109
+110 1 3 3 1 0 109 110
+111 1 3 3 1 0 110 111
+112 1 3 3 1 0 111 112
+113 1 3 3 1 0 112 113
+114 1 3 3 1 0 113 114
+115 1 3 3 1 0 114 115
+116 1 3 3 1 0 115 116
+117 1 3 3 1 0 116 117
+118 1 3 3 1 0 117 118
+119 1 3 3 1 0 118 119
+120 1 3 3 1 0 119 120
+121 1 3 3 1 0 120 121
+122 1 3 3 1 0 121 122
+123 1 3 3 1 0 122 123
+124 1 3 3 1 0 123 124
+125 1 3 3 1 0 124 125
+126 1 3 3 1 0 125 126
+127 1 3 3 1 0 126 127
+128 1 3 3 1 0 127 128
+129 1 3 3 1 0 128 129
+130 1 3 3 1 0 129 130
+131 1 3 3 1 0 130 131
+132 1 3 3 1 0 131 132
+133 1 3 3 1 0 132 133
+134 1 3 3 1 0 133 134
+135 1 3 3 1 0 134 135
+136 1 3 3 1 0 135 136
+137 1 3 3 1 0 136 137
+138 1 3 3 1 0 137 138
+139 1 3 3 1 0 138 139
+140 1 3 3 1 0 139 140
+141 1 3 3 1 0 140 141
+142 1 3 3 1 0 141 142
+143 1 3 3 1 0 142 143
+144 1 3 3 1 0 143 144
+145 1 3 3 1 0 144 145
+146 1 3 3 1 0 145 146
+147 1 3 3 1 0 146 147
+148 1 3 3 1 0 147 148
+149 1 3 3 1 0 148 149
+150 1 3 3 1 0 149 150
+151 1 3 3 1 0 150 151
+152 1 3 3 1 0 151 152
+153 1 3 3 1 0 152 153
+154 1 3 3 1 0 153 154
+155 1 3 3 1 0 154 155
+156 1 3 3 1 0 155 156
+157 1 3 3 1 0 156 157
+158 1 3 3 1 0 157 158
+159 1 3 3 1 0 158 159
+160 1 3 3 1 0 159 160
+161 1 3 3 1 0 160 161
+162 1 3 3 1 0 161 162
+163 1 3 3 1 0 162 163
+164 1 3 3 1 0 163 164
+165 1 3 3 1 0 164 165
+166 1 3 3 1 0 165 166
+167 1 3 3 1 0 166 167
+168 1 3 3 1 0 167 168
+169 1 3 3 1 0 168 169
+170 1 3 3 1 0 169 170
+171 1 3 3 1 0 170 171
+172 1 3 3 1 0 171 172
+173 1 3 3 1 0 172 173
+174 1 3 3 1 0 173 174
+175 1 3 3 1 0 174 175
+176 1 3 3 1 0 175 176
+177 1 3 3 1 0 176 177
+178 1 3 3 1 0 177 178
+179 1 3 3 1 0 178 179
+180 1 3 3 1 0 179 180
+181 1 3 3 1 0 180 181
+182 1 3 3 1 0 181 182
+183 1 3 3 1 0 182 183
+184 1 3 3 1 0 183 184
+185 1 3 3 1 0 184 185
+186 1 3 3 1 0 185 186
+187 1 3 3 1 0 186 187
+188 1 3 3 1 0 187 188
+189 1 3 3 1 0 188 189
+190 1 3 3 1 0 189 190
+191 1 3 3 1 0 190 191
+192 1 3 3 1 0 191 192
+193 1 3 3 1 0 192 193
+194 1 3 3 1 0 193 194
+195 1 3 3 1 0 194 195
+196 1 3 3 1 0 195 196
+197 1 3 3 1 0 196 197
+198 1 3 3 1 0 197 198
+199 1 3 3 1 0 198 199
+200 1 3 3 1 0 199 200
+201 1 3 3 1 0 200 2
+$EndElements
diff --git a/Solver/dgAlgorithm.cpp b/Solver/dgAlgorithm.cpp
index e93d4164964467e15076f78e9f8ee3a1b68dbd9e..a65907abdc4ca3ce80d8101638e0abb4b0f33b09 100644
--- a/Solver/dgAlgorithm.cpp
+++ b/Solver/dgAlgorithm.cpp
@@ -296,11 +296,14 @@ void dgAlgorithm::rungeKutta (const dgConservationLaw &claw, // conservation l
Unp._data->axpy(*(sol._data));
for(int j=0; j<orderRK;j++){
- if(j){
- K._data->scale(b[j]);
- K._data->axpy(*(sol._data));
- }
- // printf("iter %d sol size = %d\n",j,sol._dataSize);
+ if(j){
+ K._data->scale(b[j]);
+ K._data->axpy(*(sol._data));
+ if (limiter){
+ limiter->apply(K, eGroups, fGroups);
+ }
+ }
+
this->residual(claw,eGroups,fGroups,bGroups,K._dataProxys,resd._dataProxys);
K._data->scale(0.);
for(int k=0;k < eGroups.size();k++) {
@@ -313,14 +316,11 @@ void dgAlgorithm::rungeKutta (const dgConservationLaw &claw, // conservation l
}
}
Unp._data->axpy(*(K._data),a[j]);
- if (limiter) limiter->apply(Unp, eGroups, fGroups);
}
-
+ if (limiter) limiter->apply(Unp, eGroups, fGroups);
for (int i=0;i<sol._dataSize;i++){
- // printf("tempSol[%d] = %g\n",i,(*tempSol._data)(i));
- // memcp
+
(*sol._data)(i)=(*Unp._data)(i);
- //if (limiter) limiter->apply(sol, eGroups, fGroups);
}
}
diff --git a/Solver/dgLimiter.h b/Solver/dgLimiter.h
index d26aeccc183131f668b841a122d2a31857d3cad9..7d8cab78fa279020485a3bd572bc475f0b8fb0e3 100644
--- a/Solver/dgLimiter.h
+++ b/Solver/dgLimiter.h
@@ -14,9 +14,9 @@ public:
std::vector<dgGroupOfFaces*> &fGroup ) = 0;
};
-class SlopeLimiter : public dgLimiter{
+class dgSlopeLimiter : public dgLimiter{
public :
- SlopeLimiter () : dgLimiter () {}
+ dgSlopeLimiter () : dgLimiter () {}
virtual bool apply ( dgDofContainer &solution,
std::vector<dgGroupOfElements*> &eGroups,
std::vector<dgGroupOfFaces*> &fGroup);
diff --git a/Solver/dgSlopeLimiter.cpp b/Solver/dgSlopeLimiter.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..379172738ef20b4ff834609078e2cd922922a6c3
--- /dev/null
+++ b/Solver/dgSlopeLimiter.cpp
@@ -0,0 +1,141 @@
+#include "dgLimiter.h"
+#include "dgGroupOfElements.h"
+#include "dgSystemOfEquations.h"
+
+
+//----------------------------------------------------------------------------------
+bool dgSlopeLimiter::apply ( dgDofContainer &solution,
+ std::vector<dgGroupOfElements*> &eGroups,
+ std::vector<dgGroupOfFaces*> &fGroups)
+{
+
+ //WARNING FOR ONLY 1 GROUP OF FACES
+ dgGroupOfFaces* group = fGroups[0];
+ fullMatrix<double> &SolLeft = *(solution._dataProxys[0]);
+ fullMatrix<double> &SolRight = *(solution._dataProxys[0]);
+ int nbFields = solution.getNbFields();
+ int totNbElems = solution.getNbElements();
+
+ // first compute max and min of all fields for all stencils
+ //----------------------------------------------------------
+
+ fullMatrix<double> MIN (totNbElems,nbFields );
+ fullMatrix<double> MAX (totNbElems, nbFields);
+ MIN.setAll ( 1.e22);
+ MAX.setAll (-1.e22);
+
+ fullMatrix<double> MEANL (totNbElems, nbFields);
+ fullMatrix<double> MEANR (totNbElems, nbFields);
+ MEANL.setAll ( 0.01);
+ MEANR.setAll ( 0.01);
+
+ int iElementL, iElementR, fSize;
+ for(int iFace=0 ; iFace<group->getNbElements();iFace++)
+ {
+ iElementL = group->getElementLeftId(iFace);
+ iElementR = group->getElementRightId(iFace);
+ if (iElementR >= 0)
+ {
+ fullMatrix<double> TempL, TempR;
+ TempL.setAsProxy(SolLeft, nbFields*iElementL, nbFields );
+ TempR.setAsProxy(SolRight, nbFields*iElementR, nbFields );
+
+ fSize = TempL.size1();
+ for (int k=0; k< nbFields; ++k)
+ {
+ double AVGL = 0;
+ double AVGR = 0;
+ for (int i=0; i<fSize; ++i)
+ {
+ AVGL += TempL(i,k);
+ AVGR += TempR(i,k);
+ }
+ AVGL /= (double) fSize;
+ AVGR /= (double) fSize;
+ MIN ( iElementL , k ) = std::min ( AVGR , MIN ( iElementL , k ) );
+ MAX ( iElementL , k ) = std::max ( AVGR , MAX ( iElementL , k ) );
+ MIN ( iElementR , k ) = std::min ( AVGL , MIN ( iElementR , k ) );
+ MAX ( iElementR , k ) = std::max ( AVGL , MAX ( iElementR , k ) );
+
+ MEANR( iElementL,k ) = AVGR;
+ MEANL( iElementR,k ) = AVGL;
+
+ }
+ }
+ else{
+ fullMatrix<double> TempL ;
+ TempL.setAsProxy(SolLeft, nbFields*iElementL, nbFields );
+ for (int k=0; k<nbFields; ++k){
+ for (int i=0; i<fSize; ++i){
+ MIN ( iElementL , k ) = std::min ( TempL(i,k) , MIN ( iElementL , k ) );
+ MAX ( iElementL , k ) = std::max ( TempL(i,k) , MAX ( iElementL , k ) );
+
+
+ double AVG = 0;
+ for (int i=0; i<fSize; ++i) AVG += TempL(i,k);
+
+ }
+ }
+ }
+ }
+ // some parallel should be done here in order to compute averages on other processors
+ //----------------------------------------------------------
+
+ //...
+
+ // then limit the solution
+ //----------------------------------------------------------
+
+ for (int iElement=0 ; iElement<totNbElems; ++iElement)
+ {
+ fullMatrix<double> Temp;
+ Temp.setAsProxy(SolRight, nbFields*iElement, nbFields );
+ for (int k=0; k<nbFields; ++k)
+ {
+ double AVG = 0.;
+ double locMax = -1.e22;
+ double locMin = 1.e22;
+ double neighMax = MAX (iElement,k);
+ double neighMin = MIN (iElement,k);
+ for (int i=0; i<fSize; ++i)
+ {
+ AVG += Temp(i,k);
+ locMax = std::max (locMax, Temp (i,k));
+ locMin = std::min (locMin, Temp (i,k));
+ }
+ AVG /= (double) fSize;
+ // printf("AVG %e LocMax %e Locmin %e\n",AVG,locMax,locMin);
+
+ //SLOPE LIMITING DG
+ //-------------------
+ for (int i=0; i<fSize; ++i)Temp(i,k) -= AVG;
+
+ double slopeLimiterValue = 1.0;
+ if (locMax != AVG && locMax > neighMax) slopeLimiterValue = (neighMax-AVG) / (locMax-AVG);
+ if (locMin != AVG && locMin < neighMin) slopeLimiterValue = std::min ( slopeLimiterValue , (AVG-neighMin) / (AVG-locMin) );
+ if (AVG < neighMin) slopeLimiterValue = 0;
+ if (AVG > neighMax) slopeLimiterValue = 0;
+
+ // if (detectDISC(iElement) == 0.0) slopeLimiterValue = 1.0; // do not limit
+// if (slopeLimiterValue != 1.0) printf("limit (%e) elem =%d \n", slopeLimiterValue, iElement);
+// slopeLimiterValue=0;
+
+ for (int i=0; i<fSize; ++i) Temp(i,k) *= slopeLimiterValue;
+
+ for (int i=0; i<fSize; ++i) Temp(i,k) += AVG;
+
+ //MINMOD DG (IF CHANGE TO THIS ADD LIMITER LINE IN RUNGEKUTTANNOPERATOR.cc)
+ //--------------------------------
+ // if (detectDISC(iElement) != 0.0){
+ // const size_t nn = Temp.size1();
+ // double alpha = 0.5; //0.5 1.0 1.3 (alpha=0.5 => MUSC Schem VAN LEER)
+ // Temp(0,k) = AVG-minmod(AVG-Temp(0,k), alpha*(AVG-MEANL(iElement,k)), alpha*(MEANR(iElement,k)-AVG));
+ // Temp(nn-1,k) = AVG+minmod(Temp(nn-1,k)-AVG, alpha*(AVG-MEANL(iElement,k)), alpha*(MEANR(iElement,k)-AVG));
+ // }
+
+
+ }
+ }
+ return true;
+
+}
diff --git a/Solver/dgSystemOfEquations.cpp b/Solver/dgSystemOfEquations.cpp
index 0b2776516d4326d0926afccb0c5369619127e0e6..cb394b03402ff76cb3f877fa9f165afb41d31edd 100644
--- a/Solver/dgSystemOfEquations.cpp
+++ b/Solver/dgSystemOfEquations.cpp
@@ -47,8 +47,10 @@ methodBinding *dgSystemOfEquations::methods[]={
new methodBindingTemplate<dgSystemOfEquations,void,dgConservationLaw*>("setConservationLaw",&dgSystemOfEquations::setConservationLaw),
new methodBindingTemplate<dgSystemOfEquations,void>("setup",&dgSystemOfEquations::setup),
new methodBindingTemplate<dgSystemOfEquations,void,std::string>("exportSolution",&dgSystemOfEquations::exportSolution),
+ new methodBindingTemplate<dgSystemOfEquations,void>("limitSolution",&dgSystemOfEquations::limitSolution),
new methodBindingTemplate<dgSystemOfEquations,void,std::string>("L2Projection",&dgSystemOfEquations::L2Projection),
new methodBindingTemplate<dgSystemOfEquations,double,double>("RK44",&dgSystemOfEquations::RK44),
+ new methodBindingTemplate<dgSystemOfEquations,double,double>("RK44_limiter",&dgSystemOfEquations::RK44_limiter),
new methodBindingTemplate<dgSystemOfEquations,double,double>("multirateRK43",&dgSystemOfEquations::multirateRK43),
0};
@@ -77,11 +79,16 @@ void dgSystemOfEquations::setup(){
double dgSystemOfEquations::RK44(double dt){
- //dgLimiter *sl = new SlopeLimiter();
_algo->rungeKutta(*_claw, _elementGroups, _faceGroups, _boundaryGroups, dt, *_solution, *_rightHandSide, NULL);
return _solution->_data->norm();
}
+double dgSystemOfEquations::RK44_limiter(double dt){
+ dgLimiter *sl = new dgSlopeLimiter();
+ _algo->rungeKutta(*_claw, _elementGroups, _faceGroups, _boundaryGroups, dt, *_solution, *_rightHandSide, sl);
+ return _solution->_data->norm();
+}
+
double dgSystemOfEquations::multirateRK43(double dt){
_algo->multirateRungeKutta(*_claw, _elementGroups, _faceGroups, _boundaryGroups, dt, *_solution, *_rightHandSide);
@@ -91,6 +98,12 @@ double dgSystemOfEquations::multirateRK43(double dt){
void dgSystemOfEquations::exportSolution(std::string outputFile){
export_solution_as_is(outputFile);
}
+void dgSystemOfEquations::limitSolution(){
+ dgLimiter *sl = new dgSlopeLimiter();
+ sl->apply(*_solution,_elementGroups,_faceGroups);
+
+ delete sl;
+}
#endif // HAVE_LUA
dgSystemOfEquations::~dgSystemOfEquations(){
diff --git a/Solver/dgSystemOfEquations.h b/Solver/dgSystemOfEquations.h
index b34256a2175efcae6cda9cf70cfbd06e0192bf47..cdcb87a1eb13c96600d9d15cd741785ba275b162 100644
--- a/Solver/dgSystemOfEquations.h
+++ b/Solver/dgSystemOfEquations.h
@@ -7,6 +7,7 @@
#include "dgAlgorithm.h"
#include "dgConservationLaw.h"
#include "Gmsh.h"
+#include "dgLimiter.h"
struct dgDofContainer {
private:
@@ -56,7 +57,9 @@ public:
dgSystemOfEquations(GModel *_gm);
void setup (); // setup the groups and allocate
void exportSolution (std::string filename); // export the solution
- double RK44 (double dt);
+ void limitSolution (); // apply the limiter on the solution
+ double RK44 (double dt);
+ double RK44_limiter (double dt);
double multirateRK43 (double dt);
void L2Projection (std::string functionName); // assign the solution to a given function