I agree with your assessment and the assessment in the Mathematica StackExchange link provided by Kuba. I also agree with Alexis comment that it is an And proposition versus an Either Or proposition. For me personally, it is also an And plus a bunch of other tools proposition. I generally have more success integrating a bunch of best in class tools to create a simulation workflow than I do using an integrated tool that does everything. Much of what I do is the simulation of highly non-Newtonian fluid flow. My preference is to use Mathematica and COMSOL to study smaller scale phenomena and other CFD packages to study industrial sized models. If you can reduce the complexity of the physics for the large-scale model with Mathematica or COMSOL, you can dramatically increase your probability of a successful simulation.
Generally, I use COMSOL and Mathematica in complementary fashions. As time marches on, there is a blurring of some of the lines between the tools. For example, Mathematica is expanding their PDE capability, while COMSOL is expanding the portability capability with the COMSOL compiler. The computable document format and Wolfram Cloud did set Mathematica apart with its ability to deploy to a broad audience.
COMSOL is definitely a more mature physics modeling platform. As such, it has some convenience features that I will add in addition to yours.
- Material library
- Boundary layer meshing for CFD/non-conformal meshes/multiple reference frames
- Solving steps
- An application library of over a thousand worked examples
- Validation and verification examples
Generally, when I am approached to do new modeling project, I state that I need 3 classes of inputs to begin. Namely, I need a description of the geometry, material properties, and operating conditions. Having a material library allows one to help specify the required inputs or potentially use a proxy material if the client cannot provide the inputs.
For highly shear thinning fluids, it is necessary to create a boundary layer mesh on surfaces otherwise the shear stresses will be overestimated. It is unclear to me that this capability exists in the Wolfram language or if the FEM implementation maintains accuracy with high aspect ratio tetrahedra.
Although COMSOL touts its fully coupled multiphysics capability, it is usually much faster and more robust if you can decouple the physics into a series of solution steps. For example, if you want to study heat transfer in a mixing vessel, you could naïvely conduct a fully coupled transient simulation starting from rest. However, experience shows that you can speed up the process over 100x by first solving the flow field with the steady-state frozen rotor approximation without heat transfer, then fix the flow field and solve the transient heat equation only in the next step and get about 90 to 95% of the correct answer. COMSOL makes it relatively easy to create this sequence of study steps (toggling checkboxes and pointing to a solution to initialize from). It should be noted that most solvers can use this strategy to reduce simulation time. It is not clear to me that this capability exists in the Wolfram language in a straightforward way, but if CFD capability is in the development plan, you need the capability to break the simulation and into a series of steps otherwise simulation times will be slow.
To be an effective modeler, it is likely that you will need to create a workflow consisting of multiple steps with different solver settings. The options space is combinatorially complex. Starting from a good initial guess that you can modify to suit your purpose is much simpler than trying to discover the optimal path from the ground up. This is where an extensive application library of solving workflows can be very beneficial. As the Mathematica modeling capability matures, perhaps they could use their curating ability to create a model lookup function similar to their formula lookup function. For example, ModelLookup[non-isothermal flow] would return several types of nonisothermal flow workflows (model data objects similar to formula data objects) to choose from.
Synergy
What I like to do is integrate multiple tools together into a productive workflow. What is usually more interesting to me is not how the tools compete, but rather how I can get them to cooperate to produce something they cannot do alone.
Flight Simulators.
One of the most successful applications of simulation has been that of a flight simulator, which have been commercially available since 1929 (Link Trainer). The main purpose of a flight simulator is to train the neural nets of the prospective pilot to the real time (i.e. interactive) dynamics of multiple controls in multiple environments. We probably can agree that a flight simulator that allowed the pilot to select the aircraft, one throttle, aileron, rudder, elevator, and flap setting and returned a picture of the aircraft flying into a mountain is essentially useless. Yet that is precisely what is being asked when a client asks for the one-off simulation.
So now, we might ask is it possible to make an interactive tool for the pilots (your customers potentially down to the operator level) to explore the effects of multiple controls using Mathematica and another slow tool? Here is one possible idea.
Slider Model with Mathematica
Interactivity
When I think of an interactive slider model, Mathematicas manipulate is one of the first tools that comes to mind.
"Baking" the Physics
Since CFD simulations are generally quite slow, I am used to postprocessing the simulation in batch mode. So, let us try to make a multidimensional flipbook using MMA in another code.
Let us create a series of GIF animations varying the liquid velocity at several gas velocities while fixing the equilibrium constant in the diffusion coefficient (in Mathematica this could take a few hours as I have not optimized the process so you may want to do this overnight). When the simulation is complete, we will import the list of GIF animations.
mDic = meshfn[config -> "Co"];
basename = "varyliquid_";
filenames =
Table[basename <> IntegerString[i, 10, 3] <> "." <> "gif", {i, 1,
20}];
Table[Export[filenames[[i]],
Monitor[
Table[
modelfn[md -> mDic, k -> 0.5, dratio -> 0.5,
pel -> yy[1, 20, 0, 50, j] // N,
peg -> yy[1, 20, 0, 50, i] // N, title -> "Co-Flow"][[3]], {j,
1, 20}
],
Grid[
{{"Total progress:",
ProgressIndicator[
Dynamic[f[j, 1, 20, 1]]]}, {"j=", {Dynamic@j}}}]
],
"AnimationRepetitions" -> \[Infinity]], {i, 1, 20}];
files = FileNames["vary*.gif", NotebookDirectory[]];
gifs = Import[#] & /@ files;
It is quite easy to create a flipbook using manipulate so that we can view the effects along the liquid velocity and gas velocity dimensions.
m = Manipulate[
Show[gifs[[i, j]]], {j, 1, 20, 1}, {i, 1, Length[gifs], 1}]

Just a word of warning. When I used the SaveDefinitions->True option, MMA ballooned up to over 15GB of RAM and took a while. There probably is a more optimal approach.
Likewise, we can conduct a similar process in COMSOL to create a series of GIF animations that can be imported and displayed in manipulate.

What I like about this approach is that almost all solvers can produce a GIF animation, but few or none that I am aware of can create a flipbook along multiple dimensions. Additionally, there could be other clever tricks that one could play with Mathematica using interpolation functions on a sampled mesh, image processing, or machine learning that could be incorporated into manipulate, but I have yet to experiment with those approaches. Now, we have a pathway to distill potentially terabytes of information into a model that an operator can interact with. I find that pretty interesting.