Natural density of set of numbers not divisible by any prime in an infinite subset

The recent paper

Matomäki, Kaisa; Shao, Xuancheng, When the sieve works. II, ZBL07207214.

gives a fairly satisfactory answer to this question in the setting of arbitrary $S$. General sieve theory gives the upper bound

$$ |D(X)| \ll X \prod_{p \in S; p \leq X} (1 - \frac{1}{p})$$

(which also matches what naive probabilistic heuristics would predict) but the matching lower bound

$$ |D(X)| \gg X \prod_{p \in S; p \leq X} (1 - \frac{1}{p}) \quad (1)$$

is not always true. For instance if $S$ consists of all the primes between $X^{1/v}$ and $X$ for some fixed $v>1$ then the RHS of (1) is about $X/v$ but the LHS is instead $v^{-v+o(v)} X$ (smooth numbers are significantly rarer than naive heuristics would predict). However in this paper (resolving a previous conjecture of Granville, Koukoulopoulos and Matomaki in the reference at the end of this answer) they show (roughly speaking) that the lower bound (1) holds if one has an inequality of the form $$ \sum_{X^{1/v} \leq p \leq X^{1/u}: p \not \in S} \frac{1}{p} \geq \frac{1+\varepsilon}{u}$$ for some $v > u > 1$ that are not too large and some $\varepsilon > 0$ (there are examples that show that this condition is close to best possible); here the implied constant in (1) are allowed to depend on $u,v,\varepsilon$ and is basically of the form $v^{-e^{-1/u} v}$. Basically this condition is asserting that $S$ doesn't end up containing the majority of all primes between $X^{1/v}$ and $X^{1/u}$ for some bounded $u,v$, as this can lead to dramatic reductions in the size of $D(X)$.

In the situation where $S$ has natural density $\alpha < 1$ relative to the primes, summation by parts will give an asymptotic $\sum_{X^{1/v} \leq p \leq X^{1/u}; p \not \in S} \frac{1}{p} = (1-\alpha) \log(v/u) + o(1)$ as $X \to \infty$ keeping $u,v$ fixed, and so by choosing $u,v$ appropriately one can invoke the Matomaki-Shao theorem and obtain the matching lower bound (1) (presumably this can also be established by earlier results than the Matomaki-Shao paper). This already implies that $|D(X)| = X / \log^{\alpha+o(1)} X$ and presumably for your specific set $S$ you may be able to sharpen the $o(1)$ error term here using more quantitative versions of the Chebotarev density theorem.

For arbitrary $S$, one has a logarithmic version

$$ \sum_{n \in D(X)} \frac{1}{n} \asymp \log X \prod_{p \in S; p \leq X} (1-\frac{1}{p});$$

see Lemma 2.1 of

Granville, Andrew; Koukoulopoulos, Dimitris; Matomäki, Kaisa, When the sieve works, Duke Math. J. 164, No. 10, 1935-1969 (2015). ZBL1326.11055.

Call a set $S$ of primes "Frobenian" if there is a finite Galois extension $K/\mathbf Q$ and a union of conjugacy classes $H$ in ${\rm Gal}(K/\mathbf Q)$ such that $S$ is equal to the set of primes $p$ that are unramified in $K$ for which the Frobenius conjugacy class of $p$ in ${\rm Gal}(K/\mathbf Q)$ lies in $H$ or, more generally, $S$ is that set of primes with finitely many exceptions. An example could be the set of all primes, with finitely many exceptions, that belong to a union of arithmetic progressions modulo $m$ that are all relatively prime to $m$, such as all primes that are $1 \bmod 4$ without the primes $5$ and $29$ or the set of all primes that are $1 \bmod 4$ together with the prime $7$.

From the MO question here I was led to a paper by Serre here or here that addresses your question in Theorem 2.4(a). Suppose the Frobenian set of primes $S$ has natural density $\alpha$ where $\alpha > 0$, so it also has "regular density $\alpha$ in the sense of Delange" (Serre's defintion (1.3)). The finitely many exceptional primes allowed in the definition of Frobenian primes above corresponds to Serre's "sufficiently large" comment in his property (1.4)(c).

Your $D(X)$ is the complement of the set $E(X) = \{n \leq X : n \text{ has some prime factor } p \in S\}$. Let $E$ be the set of all positive integers with a prime factor in $S$, so $E(X) = \{n \leq X : n \in E\}$. The set $E$ is "multiplicative": a product of relatively prime positive integers is in $E$ if and only if one of the two integers is in $E$. Then $E$ fits the hypothesis of Theorem 2.4 and your $D(X)$ is what Serre calls $E'(X)$, so Theorem 2.4(a) says for $0 < \alpha < 1$ that $D(X) = E'(X) \sim cX/(\log X)^\alpha$ for some $c > 0$, and a formula for $c$ is in equation (2.6) of Serre's paper. If $\alpha = 1$ then Theorem 2.4(b) says $D(X) = E'(X) = O(X^{1-\delta})$ for some $\delta > 0$.

Under a very slight strengthening of the hypothesis, this can be done by elementary methods, namely the "Wirsing–Odoni method". The following version is Proposition 4 in this paper of mine with Finch and Sebah (I'm simplifying the hypotheses for clarity here):

Let $f$ be a multiplicative function satisfying $0\le f(n)\le 1$ for all $n$. Suppose that there exist real numbers $\xi>0$ and $0<\beta<1$ such that $$ \sum_{p<P}f(p)=\xi \frac P{\log P}+O\bigg( \frac P{(\log P)^{1+\beta}}\bigg) \tag{1} $$ as $P\rightarrow \infty$. Then the product over all primes $$ C_f=\frac1{\Gamma(\xi)} \prod_{p} \bigg( 1+\frac{f(p)}p+\frac{f(p^2)}{p^2}+\frac{f(p^3)}{p^3}+\cdots \bigg) \bigg( 1-\frac1p \bigg)^\xi $$ converges (hence is positive), and $$ \sum_{n<N}f(n)=C_fN(\log N)^{\xi -1}+O_f\big( N(\log N)^{\xi -1-\beta}\big) $$ as $N\rightarrow \infty $.

In this case, $f(p^r)$ equals $1$ if $p\notin S$ and equals $0$ if $p\in S$; then $\xi=1-\alpha$ by the assumption of the natural relative density of $S$ (we see that we need a slightly quantitative statement of that density in terms of the error term in $(1)$). In this case $$ C_f=\frac1{\Gamma(1-\alpha)} \prod_{p} \bigg( 1-\frac{1-1_S(p)}p\bigg)^{-1} \bigg( 1-\frac1p \bigg)^{1-\alpha}, $$ where $1_S(p)$ equals $1$ if $p\in S$ and equals $0$ if $p\notin S$. Notice that this is completely independent of any algebraic properties of $S$.